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SZCZECIŃSKIE ROCZNIKI NAUKOWE - ANNALES SCIENTIARUM STETINENSES 


1991 L VI, z. 1, s. 7-28 


STANISŁAW BAŃKA 


O WYZNACZANIU KANONICZNYCH REALIZACJI 
UKŁADÓW OSOBLIWYCH 
DLA NIEWŁAŚCIWYCH MACIERZY TRANSMITANCJI 
OPISANYCH WIELOMIANOWYMI 
REPREZENTACJAMI FRAKCYJNYMI* 


Instytut Automatyki Przemysłowej Politechniki Szczecińskiej, 
ul. Gen. Sikorskiego 37, 70-313 Szczecin 


Streszczenie. W pracy omówiono zagadnienia związane z określaniem kano- 
nicznych realizacji wielowymiarowych liniowych układów dynamicznych (dyskret- 
nych i ciągłych), opisywanych w dziedzinach częstotliwościowych (operatorowych) 
macierzami niewłaściwych transmitancji wymiernych. Realizacje układów wyzna- 
czane są w postaci uogólnionych (osobliwych) opisów w dziedzinie czasowej. Podano 
twierdzenia i wnioski określające warunki istnienia sterowalnych i/lub obserwo- 
walnych (minimalnych) realizacji kanonicznych, oraz podano sposób ich wyzna- 
czania, jeśli opisy układów dane są w postaci „jednostronnych” wielomianowych 
reprezentacji frakcyjnych (ułamkowych). Sposób postępowania zilustrowano przy- 
kładem obliczeniowym dokonanym dla wielowymiarowego układu ciągłego. 


Słowa kluczowe: wielowymiarowe układy dynamiczne — macierze 
wielomianowe — sterowalność i obserwowalność układów liniowych. 


WSTĘP 


W pracy rozpatruje się zagadnienia związane z wyznaczaniem realizacji 
w dziedzinie czasowej, wielowymiarowych liniowych układów dynamicznych 
(ciągłych lub dyskretnych), opisywanych macierzami niewłaściwych trans- 
mitancji w dziedzinach częstotliwościowych. Wyznaczane realizacje czasowe 
układów mają postać uwikłaną, w których lewa strona uogólnionych równań 
stanu zapisywana jest w postaci Eż (t) (lub Ez,,, — dla układów dyskretnych), 
gdzie E jest pewną (osobliwą) macierzą kwadratową. Tego rodzaju uwikłana 
forma opisów układów (niekoniecznie osobliwych) powstaje w sposób 
naturalny podczas analizy wielu procesów dynamicznych, m.in. 
w rozbudowanych układych elektrycznych i elektronicznych. W przypadku 
układów dużych, macierz E często może być osobliwa. 


* Praca była częściowo wykonywana w ramach resortowego problemu naukowo-badawczego 
RP.L.02. 


Uogólnione opisy wielowymiarowych układów liniowych można spotkać 
w wielu pracach, gdzie noszą one nazwy układów uwikłanych lub 
uogólnionych (z ang. implicit, generalized state-space, semistate lub descriptor 
systems), i jeśli E jest macierzą osobliwą — układów osobliwych (z ang. singular 
systems) [1, 2, 4-7, 11-12, 14-16, 19, 21-22, 24 i 26]. Uogólnionym opisom 
osobliwych, liniowych układów dynamicznych, odpowiadają macierze wy- 
miernych transmitancji w dziedzinach  częstotliwościowych, względem 
operatorów s=0+jw € 6 (lub z=pe'* e € dla układów dyskretnych), gdzie 
% jest ciałem liczb zespolonych. 

Podobnie jak w przypadku klasycznych (nieosobliwych) układów dyna- 
micznych, ważnym zagadnieniem jest określenie warunków dla których 
otrzymywane (uogólnione) realizacje układów osobliwych będą sterowalne 
i/lub obserwowalne w sensie impulsowym. 

W pracy dokonuje się, na podstawie najnowszej literatury światowej, 
prezentacji i częściowo uogólnień twierdzeń i wniosków umożliwiających 
wyznaczanie impulsowo sterowalnych i/lub obserwowalnych realizacji ukła- 
dów wielowymiarowych, opisywanych macierzami niewłaściwych transmitancji 
wymiernych, danych w postaci wielomianowych reprezentacji [rakcyjnych 
(ułamkowych). Wszystkie twierdzenia i wnioski są przytoczone bez dowodów, 
lecz z odnośnikami do cytowanych pozycji literatury, w których dowody te są 
zamieszczone. Dowody twierdzeń i wniosków „dualnych” pominięto ze 
względu na pełne analogie postępowania dowodowego w stosunku do 
dowodów zamieszczonych w cytowanych pracach. Sposób wyznaczania 
realizacji zilustrowano przykładem obliczeń dla wielowymiarowego układu 


ciągłego. 


WIELOMIANOWE OPISY FRAKCYJNE 
UKŁADÓW WIELOWYMIAROWYCH 
W DZIEDZINACH CZĘSTOTLIWOŚCIOWYCH 


Analizowane w pracy, wielowymiarowe układy dynamiczne opisywane 
w dziedzinach częstotliwościowych s=ó6 + jo (lub z=pe/"), za pomocą macierzy 
transmitancji wymiernych, przedstawiane są w postaci wielomianowych 
reprezentacji frakcyjnych: 


— prawostronnych 
T(s) =B, (s) Ax '(5) I,, (1) 
— oraz lewostronnych 
T(s) = 1, Az '(5) B, (5), (2) 
którym odpowiadają macierzowe równania wielomianowe w postaci: 
A, (5) xp=U 3) 


y=B,(5)x, (4) 


oraz 


4>(5)y=B;(S)u, (5) 


gdzie: ue R", ye! są wektorami sygnałów wejściowych i wyjściowych 
układu (Ściślej, ich transformat dla przebiegów nielosowych), natomiast 
x,edR*, k<n, gdzie n=deg |4A, (*)| lub n=deg|4, (*)| — jest wektorem 
sygnałów zwanych częściowym wektorem stanu lub pseudostanem układu (z 
ang. partial state. Podobne postacie mają opisy dla układów dyskretnych 
w dziedzinie częstotliwościowej z. 

Macierze transmitancji (1) i (2) są z założenia niewłaściwe, tzn. istnieje (co 
najmniej jeden) kolumnowy (wierszowy) wskaźnik j=1,2,...,m, (i=l1, 2, ..., |), 
dla którego 


deg,, B,(*)> deg.,[4,(*)] (deg,, B, (-)> deg, [A> (*)]). (6) 


Macierze wielomianowe: licznikowa B,(-)eX[:]""” i mianownikowa 
A,6)e%[:]"*""”, mogą nie być prawostronnie względnie pierwsze (p.w.p.), 
tzn. może istnieć nieunimodularny prawy dzielnik tych macierzy. Podobnie, 
macierze A;(-)eR[-]'"”" i B,(:)eR[:]'""" nie muszą być lewostronnie 
względnie pierwsze (l.w.p.), przy czym macierze mianownikowe A; (*)i A(*) są 
z założenia macierzami nieosobliwymi (zwykle, kolumnowo lub wierszowo 
właściwymi”). 

Symbolem % [:]'*" oznaczone są zbiory macierzy wielomianowych wzglę- 
dem zmiennych zespolonych s e % (lub z e ©) i współczynnikami należącymi 
do ciała liczb rzeczywistych A. Jeśli l=m>1, to R[-]'*' są niekomutatyw- 
nymi pierścieniami nad ciałem liczb rzeczywistych %, z określonymi operac- 
jami dodawania i mnożenia (nieprzemiennego), w których elementami neutral- 
nymi są: macierz zerowa 0, dla dodawania i macierz jednostkowa I, dla 
mnożenia. 

W przypadku, gdy wyznacznikami pewnych macierzy wielomianowych 
U(:)em[:]'' są dowolne liczby rzeczywiste różne od zera, to są one 
jednostkami pierścienia (elementami odwracalnymi) i nazywane macierzami 
unimodularnymi. 

Przyjęte opisy w postaci równań (1) i (2), są szczególnym przypadkiem 
ogólniejszego opisu wielomianowego 


P(:)x,=Q(")u, (7) 

y=R(')x,+V(")u, (8) 
któremu odpowiadają macierze wymiernych transmitancji określone zależ- 
nością 


T()=RG)P"'()2()+V()=T()+ FC). (©) 


* Macierz wielomianową A(')eR[']'*"" nazywamy kolumnowo (wierszowo) właściwą, 
jeżeli macierz liczbowa [,,(UT,,) utworzona ze współczynników wielomianów stojących przy 
najwyższych kolumnowych (wierszowych) stopniach tej macierzy, posiada pełny rząd r=min (I, m). 
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gdzie T,(-:)=R(-)P"'(-)Q(*) stanowi „dwustronną” frakcyjną lub sfaktoryzo- 
waną postać wymiernej części macierzy transmitancji (9), natomiast 
T,(:)=V(*) - jej część wielomianową. 

Wynika to z następującej definicji równoważności Fuhrmanna oraz związa- 
nego z nią twierdzenia, dla układów danych w postaci wielomianowych opisów 
frakcyjnych (7) - (9) [8, 9, 18]. 


Definicja 1. 


Dwie wielomianowe reprezentacje frakcyjne, macierzy wymiernych trans- 
mitancji, liniowego układu wielowymiarowego 


T()=R()P"'(-)Q()+V()=R,()Pr'()2,0)+F,0), (10) 


gdzie: R(*)P""(-)Q(*)iR,(*)P;"(:)Q,(*) są ściśle właściwymi wymiernymi 
częściami, natomiast V(:)=V,(*) — jej częściami wielomianowymi, są równo- 
ważne w sensie Fuhrmanna, jeśli istnieją macierze wielomianowe X,(*), X,(*), 
Y,(:) i Y,(*) spełniające zależność 


R Ą | P() ro Fl P,() sa FR | (1) 
X6)1I]||-RG) FC) -R, 6) V,6) | 0 I | 

przy czym X,(:)em[:]""*iP,(:)eR[-]*** są l.w.p., natomiast P(*)e 
em [-]**' i Y,(;)em[:]"'** są p.w.p. Wymiary: (kxk) macierzy P(*) 


i (k, xk,) macierzy P,(*), mogą być różne, oraz k i k, nie muszą spełniać 
warunków Rosenbrocka kż>n=deg|P(*)| i k,żn=deg|P, (*)|. 


Twierdzenie 1 [8], [9] 


Jeśli dwie reprezentacje frakcyjne, Ściśle właściwej części macierzy trans- 
mitancji (10), spełniają warunki równoważności Fuhrmanna (11) (zV(:)= 
=V,(*)=0), to ich kanoniczne (sterowalne i obserwowalne) realizacje (mini- 
malne) w dziedzinie czasowej, w przestrzeni stanów xEe%" i x,E%”, są 
podobne (tzn. wektory stanu dwóch różnych realizacji (minimalnych) związane 
są zależnością x,=Mx, gdzie M jest nieosobliwą macierzą przekształcenia 
podobieństwa), przy czym 

(i) realizacje są obserwowalne, jeśliR(-) i P(-:) oraz R,(*) i P,są p.w.p., 
(ii) realizacje są sterowalne, jeśli P(-) i Q(*) oraz P,(:) i Q,(-) są |.w.p., 
(iii) realizacje są minimalne, jeśli spełnione są warunki (i) i (ii), 

wówczas 


|P S)|=e|P, 6)|=f]s 1, — 4| (12) 
lub 

|P(z)|=2|P,(z)|=Alz7,— 4|, (13) 
gdzie « i £ są liczbami rzeczywistymi różnymi od zera. Dowód Twierdzenia 


1 (które w oryginalnej postaci było sformułowane dla układów dyskretnych, 
opisywanych właściwymi macierzami transmitancji przy V(*)=V,(*)=D), 
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przedstawiony jest w [8] i [9]. Dodatkową interpretację twierdzenia oraz 
związki z innymi definicjami równoważności układów (nieosobliwych), można 
znaleźć dodatkowo w [17], [18] i [20]. Należy jednak zauważyć, że wymierna 
część T.(*:)=R(*)P"'(*)Q(*) macierzy transmitancji (9), może zawierać w so- 
bie pewne „składowe” części wielomianowej, nie ujęte w T,(*)=V(*) (nawet 
jeśli R(-:)P"'(:)iP"*(*)Q(*) będą oddzielnie, macierzami transmitancji ściśle 
właściwymi). Stąd, w tym zapisie nie zawsze udaje się explicite oddzielić część 
ściśle właściwą macierzy transmitancji, od jej części wielomianowej [25]. 

W przypadku rozpatrywanych w pracy, „jednostronnych” opisów frakcyj- 
nych (1) i (2), problem ten nie występuje i prawostronne dzielenie macierzy 
B,(*) przez A;(*), oraz lewostronne dzielenie macierzy B,(*) przez A, (=), 
pozwalają jednoznacznie* określić część wymierną ściśle właściwą i część 
wielomianową zgodnie z zależnościami: 


T(-)=B,;6)4r'()=B,()41'€)+D()=T70)+T(), (14) 
lub 
T()=4; '(0B,(:)=42'(*)B,(:)+D()=T,()+ T(), (15) 
w których macierze B,(*), A,(*) i B,(*), A,(*) spełniają nierówności 
deg,,B,(*)<deg,,[4,(*)], J=1, 2, ..., m (16) 
lub 
deg,, B,(*)<deg,, [4„(*)], i=1, 2,..., 1 (17) 
i 
T()=D(). (18) 


UOGÓLNIONY OPIS WIELOWYMIAROWYCH 
LINIOWYCH UKŁADÓW DYNAMICZNYCH 
W DZIEDZINIE CZASOWEJ 


Niewłaściwym macierzom transmitancji w dziedzinie częstotliwościowej, 
odpowiadają opisy: 
— w formie równań stanu (dla czasu ciągłego t>0) 


x(t)=4,x(t)+B;, u(t); x,=x(to), (19) 
y()=C,x(t)+Dou()+D,u(t)+ ...+D,u"(t), (20) 
gdzie x e W”, ueR"iy e R/' są odpowiednio wektorami: stanu, wejść i wyjść 


układu, lub 


* Należy jednak pamiętać, że dzielenie przez nieregularną macierz mianownikową (tj. taką, 
która posiada nierówne stopnie kolumnowe lub wierszowe) może dawać wynik niejednoznaczny, 
o ile nie zastosuje się odpowiedniej metody dzielenia, takiej jak np. w [3]. 
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— w formie uogólnionego opisu w postaci równań uwikłanych 
Ez (t()=4z (t) + Bu(t), (21) 
y(t)=Cz(t), (22) 
gdzie, z e 8%, q>n — jest wektorem zmiennych opisujących (z ang. descriptor 
variables, semi-state lub generalized state-space variables). W przypadku ukła- 


dów dyskretnych, będą to odpowiednio zależności: 
— w postaci równań stanu (dla t=t,+n4t; n=0, 1, 2, ...) 


X,+1ZArX,PByu,; x9=x(t0), (23) 
yn=C;x,+Dou,+Dyu,,1+ --- +D,u,;, (24) 

— lub w postaci uwikłanej 
Ez,,,=Az,+Bu,, (25) 
y,=CZz,, (26) 


gdzie z e 3%, q>n — jest wektorem zmiennych opisujących dla układów 
dyskretnych. 

Jeśli w równaniach wyjść (20) lub (24) macierze D,;+0 dla i=1, 2, ..., r, to 
układy przez nie opisywane, nie spełniają postulatu przyczynowości i są 
układami osobliwymi. Wtedy odpowiadające im opisy uogólnione, posiadają 
osobliwe macierze E w równaniach (21) lub (25). 

Wśród uogólnionych układów osobliwych wyróżnia się: 
— nieregularne układy osobliwe, jeśli 


1E—A|=0, (27) 
— i regularne układy osobliwe, jeśli 
1 E—A|70, (28) 
przy czym 
deg|AE—A|=q. (29) 


Osobliwość układów objawia się m.in. tym, że przebiegi zmiennych opisują- 
cych z(t) i sygnałów wyjściowych y(t) układów ciągłych, mogą zawierać 
składowe impulsowe, mimo że sygnały u(t) nie będą zawierać żadnych 
impulsów (choć mogą zmieniać się skokowo w pewnych chwilach czasowych t). 
Inne cechy układów osobliwych, takie jak: niejednoznaczność trajektorii 
układów dyskretnych w uogólnionej przestrzeni stanów z, nieakceptowalność 
pewnych sygnałów wejściowych, antycypatywność zmiennych opisujących 
z względem sygnałów wymuszających, oraz niedopuszczalność pewnych wa- 
runków początkowych układu, opisywane są ostatnio w wielu publikacjach, 
m.in. w [16] i [2], oraz cytowanej tam literaturze. 

Podstawowym problemem rozpatrywanym w pracy, jest wyznaczanie 
realizacji niewłaściwych macierzy transmitancji w formie uogólnionych opisów, 
regularnych układów osobliwych, które dla układów ciągłych przyjmują 
następującą postać kanoniczną Weinerstrassa: [4, 6] 


110 |[z1 ©] [41 01[z, 0] [B, 
E | Ę o|-lo 5. E 4 |*|a; |:0> (30) 
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>0=1E, Cz] | (|. (31) 


gdzie z, e%" i z, e%7, oraz E, jest macierzą nilpotentną z indeksem 
kś4p=q—n, natomiast macierze A,, B,, C,, E,, B, i C,, są macierzami 
liczbowymi odpowiednio o wymiarach (nxn), (nx m), (Ixn), (px p), (p xm) 
i (Ix p). 

Przez analogię podobną postać opisu można zastosować dla układów dyskret- 
nych (25) — (26). 

Podobnie jak w przypadku układów nieosobliwych, podstawowym zagad- 
nieniem jest stwierdzenie czy wyznaczane osobliwe realizacje niewłaściwych 
macierzy transmitancji, są sterowalne i/lub obserwowalne, w określonym 
„impulsowym” sensie, oraz określenie warunków dla których impulsowo 
sterowalne i obserwowalne realizacje istnieją. 

Istnieje kilka sposobów definiowania impulsowej sterowalności i obser- 
wowalności układów osobliwych, które zależą od założonej postaci opisu 
układu, od dziedziny w której układy są analizowane, a także od przyjętych 
koncepcji ich sformułowania. Dla układów opisywanych w dziedzinach często- 
tliwościowych, definicje impulsowej sterowalności i obserwowalności wywodzi 
się od koncepcji „nieskończonych” zer odsprzęgających [19], lub poprzez 
określenie warunków dotyczących rzędów odpowiednich macierzy, w oparciu 
o macierz systemową układu, dla uogólnionego opisu w dziedzinach częstotli- 
wościowych s (lub z) [1, 21, 22]. 

W przypadku układów analizowanych w dziedzinie czasowej, ścisłe 
definicje impulsowej sterowalności i obserwowalności układów wymagają 
wprowadzenia szeregu pojęć matematycznych z teorii dystrybucji. Ponieważ 
nie są one potrzebne w dalszych częściach pracy, ograniczymy się wyłącznie do 
intuicyjnego scharakteryzowania pojęcia impulsowej sterowalności, odsyłając 
czytelnika do prac [6], [24] i [2], w których te zagadnienia są omówione 
w sposób ścisły*. 

Osobliwe układy ciągłe są impulsowo sterowalne, jeżeli przy określonych 
warunkach początkowych układu istnieją odcinkami gładkie, (nieimpulsowe) 
sygnały sterujące u(t), sprowadzające wektor zmiennych opisujących z(t) do 
zadanego zbioru dystrybucji impulsowych. Innymi słowy, jest to równoważne 
możliwości wygenerowania, za pomocą odcinkami gładkich sterowań u(t), 
maksymalnego zbioru impulsów dających możliwość usunięcia wszystkich 
impulsów z wektora zmiennych opisujących z(t), wywołanych np. zakłócenia- 
mi zewnętrznymi. Może to być osiągalne za pomocą liniowych sprzężeń 
zwrotnych, jeśli układ osobliwy jest impulsowo sterowalny. Impulsową obser- 
wowalność układów można określić w sposób dualny, z tym, że w przypadku 
nieregularnych układów osobliwych, zasada dualności nie obowiązuje. 


* Natomiast przegląd różnych definicji sterowalności i obserwowalności układów nieosob- 
liwych w dziedzinie czasowej w sensie zwykłym, można znaleźć m. in. w [13] i cytowanej tam 
literaturze. 
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REALIZACJA NIEWŁAŚCIWYCH MACIERZY TRANSMITANCJI 
W POSTACI REGULARNYCH UKŁADÓW OSOBLIWYCH 


Po dokonaniu rozdzielenia części ściśle właściwej macierzy transmitancji (1) 
lub (2), od jej części wielomianowej (poprzez prawostronne lub lewostronne 
dzielenie odpowiednich macierzy, zgodnie z zależnościami (14) lub (15)), 
problem sprowadza się do oddzielnego wyznaczenia realizacji układu dla części 
ściśle właściwej T,(*) i części wielomianowej T,()=D(-). 

Zgodnie z przyjętą dekompozycją układu (30) — (31), wyznaczane realiza- 
cje mają odpowiednio postać: 


T,(s)=C,(sI,—4,) 'B, (32) 
lub 
f,2)=C4zI,-A,)''B;, (33) 


dla części Ściśle właściwej, gdzie n=deg |A, (*)|(n=deg| A,(-)|), oraz 


T, (S)=C,(sE;—I,) 'B, (34) 
lub 
I, (2)=C, (ZE, —1,)”"B,, (35) 


dla części wielomianowej T,(-)=D(*). Przy czym p jest rzędem realizacji części 
wielomianowej, natomiast q=n+p jest rzędem realizacji czasowej całego 
osobliwego układu dynamicznego zgodnego z zapisem (30) — (31). 

Realizacja części ściśle właściwej T,(*), jest standardowym problemem 
i można go rozwiązać za pomocą jednej ze znanych metod, jak dla układów 
nieosobliwych [10, 11, 23]. Zwykle korzysta się z twierdzenia „strukturalnego” 
Wolovicha [23] — wtedy otrzymywane realizacje przyjmują II postać kanonicz- 
ną sterowalną lub obserwowalną Luenbergera-Brunowskiego — lub wykorzys- 
tuje inne postacie kanoniczne opisane np. w [10]. Natomiast realizacja części 
wielomianowej T,(*)=D(*), może być dokonana w sposób następujący. 

Macierz wielomianową (niekoniecznie kolumnowo lub wierszowo właś- 
ciwą) 


D(s)=D,+D,s+D,s*+...+D,se%[s]""" (36) 
(lub jej dyskretny odpowiednik D(z)) można zapisać odpowiednio w postaci 


D(s)=D,Sp,. (5) (37) 
lub 


D()=Sb, (8) D,, (38) 
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gdzie 
10 „0 
s0 0 
0 0 
s 0 zz 0 
0 1 0 
0 S m- 0 
Sp.(s)= (39) 
US 
00 1 
00 s 
00 Sr 
lub 
ść SAR URUJEU 0.00 0 
0.2000: Z0A07Ł 58% 87.0045 0 
5pr(s)= (40) 
000 0000 0 Is sk 


są specjalnymi macierzami wielomianowymi, natomiast D, i D, — macierzami 
liczbowymi zawierającymi odpowiednie współczynniki wielomianów macierzy 
D (s). Korzystając z tożsamości 


(I-1E)I+AE+A7E*+ ...A*" 'E' ")=I (41) 
prawdziwej dla dowolnej, nilpotentnej macierzy E o indeksie k, otrzyma się 
(SE, -I,) '=—I,-E,s—Eżs*+ ...—Ez 's'"'. (42) 


Ponieważ 
D(s)=C;,(sE,-I,) 'B,=—C,(l,+E,s+ ... +E; 's' "')B,, (43) 


łatwo zauważyć, że poszczególne wyrazy macierzowego wielomianu (36) muszą 
spełniać zależność 


D,_.,=—C,E; 'B,, dla v=l,2,...,r+l, (44) 


jeśli (34), ma być realizacją macierzy wielomianowej D (s). 
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Otrzymana zależność umożliwia wyznaczanie różnych postaci osobliwych 
realizacji dla wielomianowej części macierzy transmitancji T,(*), poprzez 
odpowiedni dobór par macierzy (E,, B,) lub (C,, E,) z nilpotentną macierzą 
E, —i stosownie do dokonanego wyboru — obliczenie pozostałej macierzy C, 
lub B,, tak aby była spełniona zależność (46) dla wszystkich wartości 
wskaźnika v. Jedną z istniejących możliwości jest taki wybór macierzy 
(£,, B,) (dualnie (Ć,, E,)), aby elementy pozostałych do wyznaczenia macierzy 
Ć,, (dualnie B,), były wprost współczynnikami macierzy wielomianowej D (5). 
Spośród wielu możliwych postaci realizacji, można wyróżnić: 

— dwie postacie impulsowo sterowalne, w których 


E,=diag[J;]; j=1,2,.-.,m, (45) 


B,=diag [b]; j=l,2,...,m (46) 


gdzie: J jj są blokowymi, kanonicznymi macierzami Jordana o wymiarach 
(k,+1) x(k,+1) odpowiednio, w postaci dolnej, lewej 
00... 0|G+O*G+u 
10. 0 
= OBO 0 (47) 
O>Brzz 31 20 
lub prawej, górnej 
10 0 |FG+0*G;+1) 
i. 0 
KG OCZ. 00 
Jd> (48) 
0 dw. 0 


którym odpowiadają 
11 4 x ([;+ 1)) 


b= (49) 


lub 
0]U=66+1) 


1 (50) 
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W przypadku pierwszym, macierz Ć, jest określona wprost jako 
Ć,=—D,, (51) 
natomiast w przypadku drugim, należy ją wyznaczyć poprzez przestawienie 


odpowiednich kolumn macierzy D,. 
Rzędem powyższych realizacji jest liczba 


p=), k,tm, (52) 
j=1 
oraz 
— dwie postacie impulsowo obserwowalne, w których 
E,=diag[J,]; i=1,2,..,1, (53) 
i 
Ć,=diag[c,]; i=1,2,..,l, (54) 


gdzie J,, są blokowymi, kanonicznymi macierzami Jordana o wymiarach 
(k,;+1)x(k,+1) w postaci prawej, górnej 


© L0ZZSBIT TWO 
0.0 1 as 0 

J,= 0 0 0 jie 0 (55) 
000 1 
0.00 0 

lub lewej, dolnej 

0 0:0 0 |kr+ Dx(ki+1) 
10 0 

J,=|0 10... 0 = 
000 0 


którym odpowiadają: 


e=|b_070 0 444, 0] 7D (57) 

lub 
G= [006 aa 1] 09, (58) 
W przypadku pierwszym, macierz B, jest określona bezpośrednio jako 
B,=—D,, (59) 


2 Szczecińskie Rocz Nauk. t VI, z 1 
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natomiast w przypadku drugim, należy ją wyznaczyć przez permutacje 
odpowiednich wierszy macierzy D,. 
Rzędem tych realizacji jest liczba 


l 
p=),k,+l. (60) 
i=1 
W powyższych zależnościach k;; j=l1,2,...,m, są kolumnowymi stop- 
niami macierzy wielomianowej D(*), natomiast k, i=1, 2, ..., l — jej stopniami 


wierszowymi. Wskażniki k, w realizacjach impulsowo sterowalnych i k, 
w realizacjach impulsowo obserwowalnych, (zwane indeksami Jordana) od- 
grywają podobną rolę jak wskaźniki sterowalności dj i wskaźniki obserwowal- 
ności d, (indeksy Kroneckera) w układach nieosobliwych [10], [23]. Ze 
względu na istniejące analogie do I i II postaci kanonicznych Luenber- 
gera-Brunowskiego układów nieosobliwych, proponuje się: aby postacie im- 
pulsowo sterowalne z prawą górną i postacie impulsowo obserwowalne z lewą 
dolną macierzą Jordana — nazywać postaciami II rodzaju. Natomiast postacie 
impulsowo sterowalne z lewą dolną i postacie impulsowo obserwowalne 
z prawą górną macierzą Jordana — postaciami I rodzaju. Należy jeszcze 
zauważyć, że przedstawione powyżej realizacje wielomianowych części macie- 
rzy transmitancji, nie są jedyne, oraz że należą one do zbioru realizacji, które 
w pracy [4] nazywają się kanonicznymi postaciami Jordana, odpowiednio, 
impulsowo sterowalnymi lub obserwowalnymi (z ang. Jordan controllable or 
observable canonical forms). 

W przypadku wyznaczania realizacji dwustronnych” wielomianowych 
opisów frakcyjnych (7) -(9), w których wydzielenie części wielomianowej 
nastręcza trudności, można wykorzystać twierdzenia i sposób wyznaczania 
realizacji (minimalnych), podane w [25]. 

Zgodnie z [1], [15], impulsowa sterowalność realizacji określonej zależnoś- 
ciami (45) — (52), jest równoważna istnieniu minora rzędu p macierzy blokowej 


[sE,—I, B,], (61) 


m 


którego wyznacznikiem jest wielomian stopnia r= ); k;, równego rzędowi 


A j=l 
macierzy E,. Dualnie, impulsowa obserwowalność realizacji (53) — (60) jest 
równoważna istnieniu minora rzędu p macierzy blokowej 


sE;—I, 
© Aj, (62) 


[i 
którego wyznacznikiem jest wielomian stopnia r= ))k,, równego rzędowi 
i=1 
macierzy E,. Realizacja impulsowo sterowalna nie musi być realizacją impul- 
sowo obserwowalną i odwrotnie, realizacja impulsowo obserwowalna może nie 
być impulsowo sterowalną. Realizacje impulsowo sterowalne i obserwowal- 
ne,są realizacjami minimalnymi. 


IE) 


Użytecznymi kryteriami impulsowej sterowalności i/lub impulsowej obser- 
wowalności (wykorzystywanymi bezpośrednio w definicjach tych pojęć) dla 
dowolnych (niekoniecznie kanonicznych) realizacji (34) lub (35),wielomianowej 
części T, (*)=D(*), jest to, aby blokowe macierze liczbowe 


[E, B,]. (63) 


E 
Hi (64) 


były macierzami pełnego rzędu równego p [1]. 

Wyznaczane realizacje w formie osobliwych opisów w uogólnionej prze- 
strzeni stanów z, obejmujące część Ściśle właściwą T,(*) i wielomianową L(): 
niewłaściwych macierzy transmitancji, mają postać (impulsowo) sterowalną 
rzędu 


i/lub 


q=n+p= 2 dj+ Ż k,+m, (65) 
ji gel 
dla (14), lub postać (impulsowo) obserwowalną rzędu 
l l 
q=n+p= ),d,+ Q,k,+l, (66) 
i=l  i=l 


dla macierzy transmitancji (15). 

Warunki, dla których otrzymywane realizacje niewłaściwych macierzy 
transmitancji danych w postaciach frakcyjnych (14) i (15) (a także (1) i (2)) są 
minimalne, można przedstawić w postaci następujących twierdzeń, praw- 
dziwych dla regularnych układów osobliwych. 


Twierdzenie 2 [7] 


Uogólniona realizacja w postaci sterowalnej rzędu q określonego zależnoś- 
cią (65), jest realizacją minimalną (impulsowo sterowalną i impulsowo obser- 
wowalną), wtedy i tylko wtedy, jeśli prawostronna wielomianowa reprezentacja 
frakcyjna macierzy transmitancji (14) jest nieredukowalna, tzn. jeśli macierze 
wielomianowe B,(*) i A,(*) (równoważnie B,(*) i 4,(*)) są p.w.p., oraz 
macierz kolumnowych współczynników wiodących T'y, macierzy D(*), posia- 
da rząd m. 


Twierdzenie 3 (dualne) 


Uogólniona realizacja w postaci obserwowalnej rzędu q określonego 
zależnością (66), jest realizacją minimalną (impulsowo obserwowalną i impul- 
sowo sterowalną), wtedy i tylko wtedy, jeśli lewostronna wielomianowa 
reprezentacja frakcyjna macierzy transmitancji (15) jest nieredukowalna, tzn. 
jeśli macierze wielomianowe A,(*) i B,(*) (ekwiwalentnie A,(*) i B,(;) są 
lw.p., oraz macierz wierszowych współczynników wiodących [„», macierzy 
D(-) posiada rząd l. 


2* 
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Twierdzenie 2 zostało udowodnione w [7], natomiast prawdziwość Twier- 
dzenia 3 wynika z zasady dualności pojęć sterowalności i obserwowalności, 
obowiązującej dla regularnych układów osobliwych [6]. Ponadto można 
wykazać słuszność twierdzeń pomocniczych, z których jedno z nich ([7], temat 
2.1.) zacytujemy w, postaci dwóch wniosków dualnych. 


Wniosek 1 [7] 


Jeśli istnieje (impulsowo) sterowalna i obserwowalna realizacja (niewłaś- 
ciwej) macierzy transmitancji (14), rzędu q określonego wzorem (65), to 
wszystkie realizacje tego samego rzędu będą również (impulsowo) sterowalne 
i obserwowalne. 


Wniosek 2 (dualny) 


Jeśli istnieje (impulsowo) obserwowalna i sterowalna realizacja (niewłaś- 
ciwej) macierzy transmitancji (15), rzędu q określonego wzorem (66), to 
wszystkie realizacje tego samego rzędu będą również (impulsowo) obserwowal- 
ne i sterowalne. 

Przedstawione powyżej rozważania i twierdzenia ilustruje następujący 
przykład obliczeniowy. 


Przykład 


Należy wyznaczyć kanoniczne (sterowalną i obserwowalną) realizacje 
wielowymiarowego liniowego układu dynamicznego o m=3 wejściach i I=2 
wyjściach, opisanego macierzą transmitancji 


s++3s84+2s*+s 1 2 
s+1 S s+1 
T()= ; 
s+1 —2_ s*+2s+1 
s s+2 s+2 


Po sprowadzeniu powyższej macierzy do postaci ułamkowej, w której diagona- 
lne macierze mianownikowe 4;,(s) i 4A,(s) utworzone są z najmniejszych 
wspólnych mianowników transmitancji, należących do odpowiednich kolumn 
i wierszy macierzy T(s), otrzyma się 


T()=B, ()Ar'()=4; '(5) B2(), 


gdzie 
s+s 0 0 
A,()=| 0  s*+2s 0 
0 0  s*+3s+2 


B, (s s” 4+3s*+28'+s* s+2 25+4 
S)= 
a s+s*+s+1 2s  s3+3s2+3s+1 
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oraz 


s +s 0 
GW] go ai 


B.(5 s” +3s*+2s3+sż s+1 25 
3 s*+28+s+2 25. „gl E2740 | 


Widoczne jest, że macierz T(s) jest macierzą transmitancji niewłaściwych 
i składać się będzie z części ściśle właściwej I, (s), oraz części wielomianowej 
T,(s). Dokonując prawostronnego dzielenia macierzy B,(s) przez A, (s) otrzy- 
ma się 


B, (3)A1*()=B,()41*(+D(9)=1,(5)+ T,(5) , 
gdzie 


20=| s+2 p 


s+l 25 s+1 


D() s +2s+1 0 0 
s)= | 
s” 0-5 
Podobnie, lewostronne dzielenie macierzy B,(s) przez A, (s) daje wynik 
4; 6)B,()=47 (5) B, (9+D()= 7,(9+ 7,(s), 


gdzie 


B,(s)= 
6) a 2s s 


z identyczną jak poprzednio macierzą wielomianową D(s). Stąd, ściśle 
właściwa część macierzy transmitancji, odpowiadająca nieosobliwej części 
realizowanego układu dynamicznego, jest określona jako 


—s s+l | 


1. 1 2 
: ; = " ż s+l S s+1 
T,()=B,(8)4r'(55=4z'(5)B,(s)= | ZAOO | |[ 


S s+2 s+2 
natomiast część wielomianowa, opisująca osobliwą część układu, jest określona 
przez macierz wielomianową 


T,(5)=D(s)=D,+D,s+D,s*+D;,s*= 


[too oj, [o © 0] [2 9 9]2,[1 0 0]. 
"66H Fóedojh Js Gl LÓD 2 
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a) Postać sterowalna 


Łatwo sprawdzić, że macierze B, (s) i A, (s) nie są p.w.p., przy czym 4, (s) 


100 
jest macierzą kolumnowo właściwą o F,=|0 1 0|, z kolumnowymi 
00 1 


stopniami (indeksami sterowalności Kroneckera) d, =2, d,=2 i dy =3 — speł- 
niającymi warunki 
deg,, [B, (s)] <deg,, A, (s)] dla j=1, 2,3. 
Na podstawie twierdzenia „strukturalnego” Wolovicha [23] i zgodnie ze 
sposobem postępowania podanym w [10], jako realizację ściśle właściwej 
części macierzy transmitancji ł,(s), otrzyma się nieosobliwy układ Ściśle 
przyczynowy, w Il-giej postaci kanonicznej sterowalnej Luenbergera-Brunow- 
skiego 
4,0=4,2,()+8, ut) 
y, (0)=Ć, z, (1) 


gdzie: 
ż,€ R"; n=2,+1,+4,=degiA, (s|=6, 
0 10 0 0 0 000 
0-10 0 0 0 100 
A 0 00 1 0 0 z 000 
A= , B,= 
0 00 -2 0 0 © 1-0 
0 00 0 0 1 00 0 
0 00 0 -2 -3 00 1 


0-1 2 14 1 
cf 1 0 2 1 Ak 
Natomiast realizację części wielomianowej T (s)=D(s), w postaci impul- 


sowo sterowalnej, otrzyma się w sposób następujący. Ponieważ kolumnowe 
10 

00 1 
ik;=1, to zgodnie z (37) i (40) można ją przedstawić w postaci 


D(s)=D,Sp.(5) 


stopnie macierzy D (s) o Ip, = , wynoszą odpowiednio k,=3,k,=0 


gdzie 
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100 

s00 

sż00 
S„GS= |" 0 01, 

0.10 

0.0 1 

0 0 s 


skąd na podstawie (45) —(52) otrzyma się opis części osobliwej układu 
w I postaci impulsowo sterowalnej 


28 ż, (t) =ź, (t) +B, u(t) 


Y2 (0) =Ć6; Źa(t), 
gdzie 


ż, e€ RP , p=k,+k,+ką+m=1 $ 


0000000 100 
1000000 000 
0100000 000 
E,=|00 1000 0|, B8=|00 0 
0000000 010 
0000000 00 1 
0000010 000 
i 
: = +0 2010405000 
C=-B--|o 91000 1] 


przy czym rząd nilpotentnej macierzy E, wynosi r=k, +k, +k;=4 (przy jej 
indeksie k=3) 


Łatwo sprawdzić, że dla otrzymanych macierzy £,, B, i Ć, spełniona jest 
równość (46), zatem otrzymany opis osobliwej części układu jest realizacją 
wielomianowej części macierzy transmitancji T,(s)=D (s) i jest to realizacja 
impulsowo sterowalna (co można potwierdzić za pomocą kryterium (63)). 
Otrzymana realizacja w postaci sterowalnej nie jest jednak impulsowo obser- 
wowalna, ponieważ zgodnie z Tw. 2, rząd macierzy I'„, nie równa się m (=3), 
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chociaż macierz D (s) jest kolumnowo właściwa (potwierdza to także niespeł- 
nione kryterium (64)). Podobnie, otrzymana w postaci sterowalnej, realizacja 
części Ściśle właściwej również nie jest obserwowalna w sensie zwykłym, 
ponieważ macierze wielomianowe B, (s) i A,(S) nie są p.w.p. 

Całkowity wymiar wektora zmiennych opisujących z e 2%, otrzymanego 
układu osobliwego, zgodnego z zapisem (30) — (31) wynosi q=n+p=6+7=13 
i w wyznaczonej postaci, realizacja ta jest impulsowo sterowalna. Jednak ze 
względu na nieobserwowalne części, osobliwą i nieosobliwą układu, nie jest to 
realizacja minimalna. 


b) Postać obserwowalna 

Łatwo wykazać, że macierze A, (s) i B, (s) są l.w.p., przy czym macierz A, (s) 
10 
o 1]> z wierszowymi stopniami (indeksami 
obserwowalności Kroneckera) d;,=2 i d,=2, spełniającymi warunki 


deg,, [B, (s)] <deg,, [A> (s)] dla i=1, 2. 


jest wierszowo właściwa o I,,= 


Podobnie jak poprzednio, na podstawie twierdzenia „strukturalnego” Wolovi- 
cha otrzyma się dualną wersję realizacji części ściśle właściwej T, (s) w II postaci 
kanonicznej obserwowalnej Luenbergera-Brunowskiego 


ż,(0)=4, z, (16)+B, u(t) 
y1(6)=C, 2, (1) 
gdzie: 
z, e%'; n=d,+d,=deg|4,(s)|=4, 


0 00 0 010 
> 1-10 0 A -1 1 2 
A,= A B,= 

0 00 0 200 

0 0 1 -2 21 


g 0100 

€,-|o 00 | 
Realizację w postaci impulsowo obserwowalnej części wielomianowej 
T,(s)=D(s) otrzyma się następująco. Zgodnie z (38) i (40), macierz wielo- 


10 


100 
stopniach k,=3 i k,=2, można przedstawić w postaci 


mianową D(s) (wierszowo niewłaściwą) o I„,= | oraz wierszowych 
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D(s)=Sp,(s5) D, 
gdzie 


Gi 
Il 
- O © FE NO 


skąd, na podstawie (53) — (60) otrzyma się opis części osobliwej układu 
w I postaci impulsowo obserwowalnej opisanej równaniami: 


E,ż,()=ż, (1)+ B,u(t), 
J2(1)=Ć, zz (1) 


gdzie 

z, EM”; p=k,+k,+1=3+2+2=7, 
0100000 1 6.0 
0010000 000 
0-.Goó [6100 2-60 

E,=|00 000 0 0|, BEHE=-B=-|100 
0000010 000 
0000 0 0 1 6001 
0000000 1-0: © 


10000 
ZACZ E 


00 
o 6] 
przy czym rząd nilpotentnej macierzy E, wynosi obecnie r=k, +k,=5, przy jej 
indeksie k=3. 

Łatwo sprawdzić, że dla otrzymanych macierzy E,, B, i Ć, spełniona jest 
równość (46) i wyznaczona w tej postaci realizacja części wielomianowej D (s) 
jest impulsowo obserwowalna. Nie jest to jednak realizacja impulsowo 
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sterowalna, ponieważ zgodnie z Tw. 3, rząd macierzy I'p, nie równa się /(=2) 
(macierz D(s) nie jest wierszowo właściwa). Potwierdza to również kry- 
terium (63). 

Rząd otrzymanego regularnego układu osobliwego w postaci obserwowal- 
nej, zgodnego z zapisem (30) — (31), wynosi obecnie q=n+p=4+7=11 
i wyznaczona realizacja jest impulsowo obserwowalna. Jednak realizacja 
ta nie jest minimalna, ponieważ posiada impulsowo niesterowalną część 
osobliwą. 

Należy zauważyć, że usunięcie największego wspólnego prawego dzielnika 
z macierzy B, (s) i A; (s) (równoważnie, B, (s) i A,(s)) przed wyznaczeniem ich 
realizacji, gwarantuje (zgodnie z pierwszą częścią wymagań Tw. 2), że nieosob- 
liwa część układu będzie jednocześnie, sterowalna i obserwowalna w sensie 
zwykłym. Tego rodzaju operacji nie można jednak przeprowadzić dla części 
wielomianowej, ponieważ nie można dokonać zmian rzędów macierzy współ- 
czynników wiodących I'p. i !p,. W rozpatrywanym przykładzie impulsowa 
nieobserwowalność realizacji w postaci sterowalnej wynika z tego, że jest za 
mało wyjść, natomiast impulsowa niesterowalność realizacji w postaci obser- 
wowalnej z tego, że macierz D(s) nie jest wierszowo właściwa. 


UWAGI KOŃCOWE 


Przedstawione w pracy rozważania umożliwiają jednoznacznie określenie 
warunków dla których wyznaczane realizacje niewłaściwych macierzy trans- 
mitancji w formie uogólnionych układów osobliwych są (impulsowo) sterowal- 
ne i/lub obserwowalne. Określają je dokładnie Twierdzenia 2 i 3 wraz 
z wnioskami. 

Sposób rozdziału zadanej macierzy (niewłaściwych) transmitancji na część 
ściśle właściwą i część wielomianową, szczególnie łatwo realizowany dla 
„jJednostronnych” wielomianowych reprezentacji frakcyjnych, umożliwia efek- 
tywne wyznaczanie kanonicznych (i minimalnych, jeśli istnieją) realizacji 
w dziedzinach czasowych, dla dowolnych macierzy transmitancji wymiernych. 
Dość istotną, numeryczną cechą przedstawionych realizacjj w postaciach 
impulsowo sterowalnych i obserwowalnych jest to, że ich parametry zawierają 
wprost wartości współczynników realizowanych macierzy wielomianowych. 
Dotyczy to również realizacji nieosobliwych części układów, z tym że po 
usunięciu odpowiednich wspólnych dzielników z macierzy licznikowych i mia- 
nownikowych zadanych macierzy transmitancji, ich macierze mianownikowe 
stają się zwykle niediagonalne. Wtedy do wyznaczenia realizacji w Il-gich 
postaciach kanonicznych Luenbergera-Brunowskiego, należy dokonać pełnego 
przebiegu obliczeń, zgodnie z zależnościami wynikającymi z twierdzenia 
„strukturalnego” Wolovicha. 


Autor pragnie podziękować recenzentom za wnikliwe uwagi, które po- 
zwoliły usunąć nieścisłości i wyjaśnić niektóre wątpliwości w czasie pisania 
i niniejszej pracy. 
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STANISŁAW BAŃKA 


ON DETERMINATION OF CANONICAL REALIZATIONS 
OF SINGULAR SYSTEMS 
FOR IMPROPER TRANSFER MATRICES GIVEN 
IN POLYNOMIAL MATRIX FRACTION DESCRIPTION 


Summary 


Some problems of impulse controllability and impulse observability of singular multivariable 
linear dynamic systems (both continuous and discrete) are considered in the paper. Analyzed 
systems are described by generalized descriptor form of state and output equations in time domain 
and by improper rational transfer matrices in frequency (s or z) domain, respectively. 

The conditions to ensure that impulse controllable and/or observable realizations of these 
transfer matrices exist are presented. Then the procedure for realization of that systems in time 
domain is given. A numerical example performed for multivariable continuous system that 
illustrates a way of employing of the results obtained is included. 


CTAHHCJIAB 5EAHbBKA 


OB OIIPEHEJJEHAM KAHOHAUUECKIHX PEAJIM3ALNMA 
OCOBEHHBIX CHCTEM 

JIA HECOBCTBEHHPBPIX MATPKN TPAHCMKTAHNUAK, 

OIIMCAHHPIX MHOTOUJEHHBIMYH OPAKNMOHHBIMM 
IIPEICTABJIEHMAMM 


Pe3rome 


B paóoTre oÓócyx/1eHbI BOIPOCHI, KACAIOLIKECH (UMNYJILCHOIi) yNpaBJIAeMOCTH M Ha- 
610 qaeMOCTH OCOÓEHHBIX MHOTOMEPDHBIX JIMACIiHBIX NAHAMUUECKUHX CHCTEM (AMHCKPETHBIX H He- 
npepBIBHBIX), OMACHIBAEMBIX B 4ACTOTHBIX (OIEpaTUBABIX) OÓJ1ACcTAx MaTpHIIAMH HECCOÓCTBEHHBIX 
PaLNioHaJIbHbIX TpARHCMHTAHIMii, A Take pEAJIHZOBAHHBIX B BHNe OGOÓLIEHHKIX YpaBAHeHNii 
COCTOAHHA H BBIXOJ|OB B OOJIACTH BpeMEHH. IIpeqCTABJIeHBI POPMYJNI H BBIBO,4BI, ONPEZIEJAŁOLIKE 
YCJIOBHA CYINECTBOBAHHA YyHIpABJIAEMLIX H(HJIM) HaOJIIONaeMBIX (MHHHMAJIEHbIX) KAHOHH4ECKHX 
pea.ruzanii, a Tae CILOCOĆ HX ONPENEJIEHKA, ECJIH ONUCAHHA CHCTEM NAHBI B BAJĄE MHOTOWIEHHBIX 
bpakuMOHABIX HpeĄCTABJIEHHi. CIOCO6 NoBEEHHA HJLIŁOCTPUDOBAH HPHUMEPOM pACYETOB, IPOH3- 
BEJIeHHBIX /UI4 MHOTOMEpHOIii HeripepBIBAOii CHCTEMBI. 
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PROBLEMATYKA WOJEWÓDZTWA SZCZECIŃSKIEGO 
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OCHRONY ŚRODOWISKA PRZYRODNICZEGO 
DO 2010 R. 


Instytut Inżynierii Wodnej Politechniki Szczecińskiej, 
al. Piastów 50, 70-311 Szczecin 


Streszczenie. Opracowano analizę zbieżności oceny stanu środowiska przy- 
rodniczego i zaopatrzenia w wodę województwie szczecińskim, przedstawioną 
w programie narodowym z rozpoznaniem wykonanym w województwie. Porównano 
koncepcję ochrony środowiska przyrodniczego i gospodarki wodnej zaproponowaną 
w programie narodowym, do głównych kierunków działania w tym zakresie, 
wypracowanych przez zespoły specjalistów z województwa szczecińskiego. 


Słowa kluczowe: ochrona atmosfery, wód powierzchniowych, podziemnych 
i gleby przed zanieczyszczeniem — gospodarka wodą, ocena stanu istniejącego 
i prognoza, koncepcja rozwiązania. 


WSTĘP 


Autorzy projektu narodowego programu ochrony środowiska przyrod- 
niczego [1] przewidują, że po sporządzeniu programu narodowego opracowa- 
ne zostaną na jego podstawie programy wojewódzkie. Powinny one precyzo- 
wać główne cele ekologiczne i sposoby ich realizacji zgodnie ze specyfiką 
przestrzenną i funkcjami społeczno-gospodarczymi województwa, a także 
zgodnie z celami strategicznymi i polityką realizacyjną, określonymi w pro- 
gramie narodowym. 

Aby narodowy program mógł być rzeczywistą podstawą programów 
wojewódzkich powinien spełniać przede wszystkim dwa z wymienionych w nim 
[1] warunków. Powinien być programem przestrzennym, uwzględniającym 
zbieżność kierunków ochrony środowiska i gospodarki wodnej w skali całego 
kraju, z działaniami w skali konkretnego regionu. Drugim z podstawowych 
warunków jest otwartość programu na okres 20 lat, oznaczająca możliwość 
jego modyfikacji co 5 lat, w oparciu o kolejne oceny stanu środowiska i wyniki 
rozwoju społeczno-gospodarczego. Celem niniejszego opracowania jest analiza 
przedstawionego projektu narodowego programu ochrony środowiska pod 
kątem jego przydatności jako podstawy do opracowania programu ochrony 
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środowiska i gospodarki wodnej województwa szczecińskiego, uwzględniające- 
go główne zadania rozwoju społeczno — gospodarczego województwa. Analiza 
obejmuje problematykę ochrony środowiska przyrodniczego z jej głównymi 
składowymi, ochroną atmosfery, wód, gleby, ochrony przed powodzią oraz 
podstawowym dla województwa szczecińskiego problemem zaopatrzenia 
w wodę ludności 'i przemysłu. 


OCENA STANU ISTNIEJĄCEGO 


Przedstawiona w programie narodowym opinia o niezahamowanej jeszcze 
tendencji zwiększania się zanieczyszczenia i degradacji środowiska przyrod- 
niczego jest w pełni adekwatna do sytuacji, jaka występuje w województwie 
szczecińskim. Ocena efektywności inwestycji, zrealizowanych w województwie 
dla ochrony środowiska, dokonana w 1987 r. przez Komisję Ochrony 
Środowiska, Gospodarki Wodnej i Leśnictwa WRN wykazała, że pomimo 
dużego zaangażowania środków w szeroki front budowy oczyszczalni ścieków 
komunalnych i przemysłowych, oraz urządzeń do oczyszczania gazów, w dal- 
szym ciągu wzrasta zanieczyszczenie środowiska przyrodniczego. 

Wymieniony stan środowiska stał się uzasadnioną podstawą zakwa- 
lifikowania terenów aglomeracji szczecińskiej i miejscowości pasa nadmor- 
skiego województwa szczecińskiego do 27 obszarów, ekologicznie zagrożonych 
na terenie kraju. Pod pojęciem terenów ekologicznie zagrożonych rozumie się 
te, na których występują potencjalne zagrożenia dla zdrowia ludzi i rozwoju 
społeczno-gospodarczego. Na obszarach tych wyodrębniono ponadto tereny 
o znacznie zdegradowanym środowisku, dla których przewiduje się podej- 
mowanie działań ochronnych z maksymalnymi preferencjami [1]. Do takich 
terenów zakwalifikowane zostały w programie narodowym tereny przyległe do 
Zakładów Chemicznych Police. 

W związku z powyższym, zagrożenie lasów w województwie szczecińskim, 
które zakwalifikowano w programie narodowym jako „Średnie” jest oceną zbyt 
optymistyczną i nieadekwatną do stanu degradacji środowiska. 42 zakłady 
przemysłowe województwa szczecińskiego, zaliczone do szczególnie uciąż- 
liwych dla powietrza atmosferycznego, emitowały w 1986 r. około 41,7 tys. ton 
pyłów i 190 tys. ton gazów, w tym 127,6 tys. ton dwutlenku siarki, 35,7 tys. ton 
tlenków azotu i 19,5 tys. ton tlenków węgla [3]. Tak znaczna emisja dwutlenku 
siarki wraz z transgranicznym napływem z terenów NRD (Schwedt, Neunbran- 
denburg, Rostock) spowodowały zniszczenie lasów na powierzchni 988 ha oraz 
zagrożenie lasów na powierzchni 2382 ha w nadleśnictwach Trzebież, Gole- 
niów i Gryfino [3]. 

Problem zanieczyszczenia powietrza w pasie nadmorskim występuje we 
wszystkich piętnastu miejscowościach-uzdrowiskach nadmorskich województ- 
wa szczecińskiego (Świnoujście, Międzyzdroje, Wisełka, Międzywodzie, Dziw- 
nów, Dziwnówek, Łukęcin, Pobierowo, Pustkowo, Trzęsacz, Rewal, Śliwno, 
Niechorze, Pogorzelica, Mrzeżyno) a nie tylko w czterech, wymienionych 
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w programie [1]. W wymienionych wyżej miejscowościach przebywa w sezonie 
około 5,5 mln osób z całego kraju [2], a wielkość ta ma tendencję rosnącą. Nie 
mniej istotnym jest stopień zanieczyszczenia powietrza na terenie aglomeracji 
szczecińskiej, gdzie w okresie zimowym normy dopuszczalnych stężeń dwutlen- 
ku siarki są przekroczone pięciokrotnie [3]. Szczególnie niebezpieczna dla 
zdrowia części mieszkańców aglomeracji jest kumulacja zanieczyszczeń pyło- 
wych i gazowych w utrzymujących się mgłach inwersyjnych, w dolinie rzeki 
Odry. 

Stan czystości wód powierzchniowych województwa szczecińskiego, przed- 
stawiony został w programie narodowym [1] w sposób niekompletny. Na 
mapie „Ochrona jakości wód”, która między innymi ilustruje stan czystości 
wód powierzchniowych w kraju, pominięty został stan czystości przybrzeżnych 
wód Bałtyku oraz wód Zalewu Szczecińskiego. Zgodnie z obowiązującymi 
przepisami, ochrona przed zanieczyszczeniem morskich wód wewnętrznych 
i przybrzeżnych znajduje się w zakresie obowiązków odpowiednich wydziałów 
urzędów wojewódzkich a nie urzędów morskich, które zajmują się regulacją 
i ochroną brzegów tych akwenów. Stąd, ocena stanu czystości i program 
ochrony tych wód przed zanieczyszczeniem powinien stanowić część składową 
programu narodowego i w większym oraz bardziej szczegółowym zakresie, 
programu wojewódzkiego. Stan czystości przybrzeżnych wód Bałtyku jest ze 
zrozumiałych powodów przedmiotem żywego zainteresowania ogółu społe- 
czeństwa, ze względu na ich funkcję rekreacyjno-leczniczą. To żywe zaintereso- 
wanie stanem czystości przybrzeżnych wód Bałtyku jest jednym z głównych 
czynników intensyfikowania działalności społecznej i edukacyjnej w zakresie 
ochrony środowiska. 

Jakość wód Zatoki Pomorskiej i Zalewu Szczecińskiego jako wód granicz- 
nych jest kontrolowana w ramach współpracy przez odpowiednie jednostki 
administracyjne strony polskiej i niemieckiej. Z badań wód morskich w pasie 
wybrzeża szczecińskiego wynika, że odpowiadają one drugiej klasie jakości 
oraz na niewielkich odcinkach trzeciej klasie (Mrzeżyno, Rewal) [3]. Stan 
czystości wód Zalewu Szczecińskiego znajduje się na pograniczu drugiej 
i trzeciej klasy. Wyniki badań wód Zalewu z okresu 1950—1986, wykazują 
tendencję wzrostu fosforanów i obfitych zakwitów sinic w powierzchniowej 
warstwie wody [4]. Wzrost odczynu wody na obszarach zakwitów powoduje, 
że staje się ona bezklasowa. 

Stan czystości rzeki Odry na odcinku aglomeracji szczecińskiej został 
oceniony w programie narodowym zbyt optymistycznie. W wyniku dopływu 
do rzeki dużych ilości ścieków miejskich i przemysłowych, czystość wody nie 
odpowiada 2-giej klasie czystości, jak to oznaczono na mapie „Ochrona jakości 
wód” [1]. Badania wody wykonane w okresie 1976-1984 wykazują stały 
wzrost zanieczyszczenień w rzece Odrze [3]. Stężenia manganu, amoniaku 
i tlenu rozpuszczonego odpowiadają trzeciej klasie czystości, a barwa wody 
i stężenia fosforanów przekraczają dopuszczalne wielkości trzeciej klasy. 

Przedstawiona na mapie [1] ocena ilości i jakości ścieków komunalnych 
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i przemysłowych odprowadzanych z miast Szczecin-Police do Odry jest 
korzystniejsza od stanu rzeczywistego. Wielkość odpływu ścieków wymagają- 
cych oczyszczania jest około dwa razy większa (—200 mln m? rocznie) od 
przedstawionej na mapie [3]. Ilość ścieków przemysłowych i komunalnych, 
oczyszczanych mechanicznie stanowi zaledwie około 20% całości, co nie jest 
zgodne z przedstawioną na mapie [1] wielkością około 60% ścieków oczysz- 
czanych mechanicznie. 

W ocenie zanieczyszczenia środowiska odpadami przemysłowymi i gos- 
podarczo-bytowymi, przedstawionej w programie narodowym zwrócono uwa- 
gę na ilościowe nagromadzenie 71% odpadów w trzech województwach 
(katowickie, legnickie i wałbrzyskie). Istotnym czynnikiem zagrożenia środowi- 
ska jest poza charakterystyką ilościową, stopień szkodliwości odpadów, 
opanowanie technologii unieszkodliwiania i zagospodarowania odpadów, 
a także możliwość pozyskania terenów o odpowiednich warunkach hydro- 
geologicznych do składowania odpadów. 

Do odpadów szczególnie niebezpiecznych dla środowiska należą fosfogipsy 
składowane w województwie szczecińskim w Policach w ilości około 5 mln ton 
rocznie, na składowisku o wielkości 28 mln ton [3]. Opracowanie technologii 
unieszkodliwiania i zagospodarowania fosfogipsów znajduje się stale jeszcze 
w fazie koncepcyjnej. Równie istotne zagrożenie dla środowiska stwarza 
w województwie szczecińskim gnojowica, powstająca w 110 dużych fermach 
bezściołowych. Roczna ilość gnojowicy wynosi 1,3 mln m* [3]. Jej unieszkod- 
liwianie i użytkowanie przeprowadzane jest w stopniu niedostatecznym, 
powodując wzrastającą degradację gleb użytkowanych rolniczo oraz zanieczy- 
szczenie wód powierzchniowych i gruntowych. Brak terenów odpowiednich 
pod składowiska odpadów komunalnych stałych i płynnych stwarza za- 
grożenie dla środowiska zwłaszcza na terenach aglomeracji szczecińskiej 
i wokół miejscowości pasa nadmorskiego. 

Zakres programu narodowego jest szerszy niż to wynika z tytułu opraco- 
wania. Poza ochroną atmosfery, wód, gleby, zasobów przyrody i kopalin, 
ochroną przed hałasem i promieniowaniem oraz ochroną przed powodzią, 
program obejmuje elementy gospodarki wodnej, nie wymienionej w tytule. 
Należą do niej: zaopatrzenie w wodę miast, wsi i przemysłu, zbiorniki wodne, 
przerzuty wody, nawodnienie użytków rolnych regulacja rzek, stopnie wodne 
i kanały żeglowne oraz elektrownie wodne. Ocena stanu gospodarki wodnej 
została przedstawiona w programie narodowym w sposób bardzo skrócony 
(około 1,5 strony i poglądowa mapa) i nieproporcjonalny do zakresu pro- 
blemów ochrony środowiska. Pozostaje to w sprzeczności z założeniem 
o kompleksowym i interdyscyplinarnym charakterze programu. Równorzęd- 
ność problematyki gospodarki wodnej, a w tym w szczególności zaopatrzenie 
w wodę gospodarki komunalnej, przemysłu i rolnictwa występuje szczególnie 
dobitnie w województwie szczecińskim. 

Program narodowy określa wielkość potencjalnego globalnego deficytu 
wody do roku 1986 na około 17% poboru wody oraz podaje, że w 25% miast, 
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w tym w 22 miastach wojewódzkich, występują stale lub okresowe deficyty 
wody. Brak jest w programie [1] istotnego elementu oceny stanu zaopatrzenia 
w wodę jakim jest określenie regionów i miast, gdzie występują maksymalne 
deficyty wody i gdzie brak wody stanowi istotną barierę dla prawidłowej 
gospodarki. 

Do wyżej wymienionych regionów, w znaczeniu ogólnokrajowym należy 
strefa nadmorska województwa szczecińskiego, z 14 miejscowościami-kąpieli- 
skami i miastem Świnoujście. W okresie sezonu w wymienionym regionie 
przebywa około 5,5 mln wczasowiczów i kuracjuszy z całego kraju. W ostat- 
nich dwóch latach (1987—1988) liczba ludności sezonowej wykazuje tendencję 
wzrostową w związku z występującym zanieczyszczeniem morskich wód 
przybrzeżnych w województwie gdańskim. Deficyt wody w mieście Świnoujście 
wynosi obecnie około 100% wydajności istniejącego wodociągu, a w pozos- 
tałych miejscowościach-kąpieliskach około 67% wydatku istniejących wodo- 
ciągów [3]. 

Przeprowadzone badania wykazują praktyczną niemożliwość zwiększenia 
wydajności istniejących ujęć wód gruntowych na wyspach Uznam i Wolin, 
ze względu na brak dostatecznych zasobów wody słodkiej. Niewielkie zasoby 
w rejonie Dargobądzia na wyspie Wolin nie mogą być wykorzystane do 
celów wodociągowych ze względu na ochronę Wolińskiego Parku Naro- 
dowego. Stwierdza się również występującą tendencję do stałego pogarszania 
się jakości wody w istniejących ujęciach wody gruntowej w pasie nadmorskim, 
co jest wynikiem nadmiernej eksploatacji ujęć oraz braku odpowiednich 
stref ochronnych. Ze względu na istotną funkcję, jaką spełnia pas nadmorski 
województwa szczecińskiego dla wypoczynku i rekreacji ludności całego 
kraju oraz niewykorzystane możliwości rozwoju funkcji sanatoryjnej dla 
kraju i zagranicy, brak wody stanowi w tym przypadku istotną barierę 
w racjonalnym wykorzystaniu walorów środowiska przyrodniczego tego 
regionu. 

Znaczące deficyty wody występują także aktualnie w systemie wodo- 
ciągowym Szczecin-Police, stanowiąc około 17% istniejącej zdolności eks- 
ploatacyjnej wodociągu [3]. Niezależnie od wymienionej wielkości deficytu 
wody występują częste awarie jedynego rurociągu magistrali wodociągowej, 
łączącej Szczecin z głównym ujęciem wody z jeziora Miedwie (50% wydajności 
wodociągu), co powoduje znaczne okresowe ograniczenie rozbioru wody 
w mieście. 


PROGRAM OCHRONY ŚRODOWISKA I GOSPODARKI WODNEJ 


W programie narodowym przyjęto jako główny cel do roku 2010 zahamo- 
wanie i odwrócenie niekorzystnych trendów w zanieczyszczeniu środowiska 
oraz doprowadzenie go do stanu, umożliwiającego skuteczną ochronę zdrowia 
człowieka, właściwe warunki egzystencji przyrody oraz rozwój gospodarczy 
kraju. 


3 Szczecińskie Rocz Nauk. t. VI, z 1 
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Dla osiągnięcia zamierzonych strategicznych celów przewiduje się w za- 
kresie ochrony atmosfery obniżyć emisję dwutlenku siarki o ponad 53% (z 4,3 
do 2,0 mln ton), co pozwoli zmniejszyć stężenie SO, na całym obszarze kraju 
poniżej normatywnej wielkości (0,032 mg/m”), a dla specjalnie chronionych 
poniżej połowy normatywnej wielkości. 

Skuteczna ochrona lasów w województwie szczecińskim, w tym Woliń- 
skiego Parku Narodowego oraz zachowanie wszelkich dopuszczalnych norm 
dla ludzi wymaga zmniejszenia emisji dwutlenku siarki o około 70% (ze 133 do 
33 tys. ton). Aby osiągnąć wymienioną redukcję SO, i równocześnie zmniejszyć 
emisję innych zanieczyszczeń atmosfery (tlenki azotu, tlenek węgla) na terenie 
województwa szczecińskiego trzeba zrealizować program obejmujący: budowę 
stacji odsiarczania i denitryfikacji spalin emitowanych z istniejących elektrowni 
(Dolna Odra, Pomorzany, Gdańska) i nowo budowanych zakładów energetyki 
cieplnej (prawobrzeżny Szczecin, Stargard, Police, Gryfice), budowę stacji 
odsiarczania spalin w kotłowniach przemysłowych i komunalnych o zdolności 
powyżej 2,093:1019J (5 Gcal), zastosowanie gazu ziemnego w miejsce węgla 
i koksu w gospodarstwach domowych, ogrodnictwach, obiektach wczasowych 
i sanatoryjnych w pasie nadmorskim, budowę stacji niszczących węglowodory 
emitowane z lakierni, malarni i stacji dystrybucji paliw [3]. 

Realizacja wojewódzkiego programu ochrony atmosfery jest uwarunkowa- 
na możliwością finansowania go z funduszu siarkowego, powstającego z opłat 
zakładów emitujących SO», uruchomieniem krajowej produkcji urządzeń do 
odsiarczania i usuwania tlenków azotu ze spalin oraz zorganizowaniem przez 
resorty przemysłu i budownictwa, specjalistycznego wykonawstwa robót budo- 
wlano-montażowych dla instalacji odsiarczania paliw oraz usuwania tlenków 
azotu ze spalin. 

Jako główne zadania w ochronie wód przed zanieczyszczeniem program 
narodowy wymienia: budowę około 2000 dużych i średnich oczyszczalni 
ścieków oraz około 8600 oczyszczalni ścieków na wsiach, uporządkowanie 
gospodarki wodno-ściekowej zakładów przemysłowych, miast i wsi, w tym 
przedsięwzięcia w zakresie zmniejszania zużycia wody w produkcji i przez 
ludność oraz nowelizację przepisów ochrony wód. 

Z mapy ilustrującej w programie narodowym jakość wód powierzch- 
niowych wynika, że jakość wody w rzece Odrze, na odcinku województwa 
szczecińskiego, w okresie 1985 --2010 pozostanie niezmiennie w drugiej klasie 
jakości. 

To założenie jest sprzeczne z wojewódzkim programem rozwoju [2] 
i koncepcją ochrony wody rzeki Odry przed zanieczyszczeniem. Istniejące 
ujęcie wody dla wodociągu szczecińskiego (Pomorzany, Zdroje), rozbudowy- 
wane (Mścięcino) oraz projektowane (Międzyodrze), które pokrywać będzie 
w przyszłości około 70% potrzeb aglomeracji pobierają, względnie pobierać 
będą wodę odrzańską, albo wodę podziemną infiltrowaną z koryta rzeki Odry. 
Od jakości wody w Odrze zależeć będzie zatem jakość wody w wymienionych 
ujęciach. W celu poprawy istniejącego stanu jakości wody w Odrze, gdzie 
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odprowadza się obecnie około 90% ścieków, wymagających oczyszczania, 
w programie wojewódzkim [3] przewiduje się do 2000 r. budowę części 
mechanicznej dużej oczyszczalni ścieków miejskich Pomorzany” oraz wyko- 
nanie części biologicznej oczyszczalni ścieków w Zdrojach. Do roku 2010 
przewiduje się budowę części mechanicznej dużej oczyszczalni ścieków miejs- 
kich „Skolwin” i części biologicznej oczyszczalni „Pomorzany”. Budowa 
wymienionych dużych oczyszczalni ścieków miejskich wymaga finansowania 
z centralnych funduszów celowych. Program ochrony wód rzeki Odry przed 
zanieczyszczeniem [3] przewiduje ponadto, że do 2000 roku zakończone 
zostaną budowy oczyszczalni ścieków przemysłowych zakładów, które powo- 
dują największe zanieczyszczenie wód (Zakłady Chemiczne „Police”, Zakłady 
Chemiczne „„Szczecin”, Zakłady Chemiczne „Chemitex-Wiskord”, Szczecińskie 
Zakłady Nawozów Fosforowych, Zakłady Papiernicze „Skolwin”, Przedsiębio- 
rstwo Połowów Dalekomorskich „Gryf” i inne). 

Mając na względzie niepełną efektywność ochrony wód Dolnej Odry przed 
zanieczyszczeniem bez zmniejszania ładunków zanieczyszczeń odprowadza- 
nych w górnym i średnim biegu rzeki, przewiduje się zorganizowanie Spółki 
Wodnej „Odra”, obejmującej głównych użytkowników rzeki, jako odbiornika 
ścieków. Wymienione przedsięwzięcia, wraz z realizacją koncepcji ograniczenia 
zrzutu zasolonych wód kopalnianych do rzeki Odry przedstawioną w pro- 
gramie narodowym, pozwalają przypuszczać, że do roku 2010 nastąpi nie tylko 
zahamowanie tendencji stałego pogarszania się jakości wody w Dolnej Odrze, 
ale poprawa jakości w stosunku do stanu istniejącego. Osiągnięcie tego celu 
jest konieczne nie tylko ze względu na potrzeby wodociągowe nowego 
budownictwa, ale również ze względu na potrzeby wodne specjalnej strefy 
ekonomicznej w Szczecinie [2]. 

Program narodowy nie ustosunkowuje się do ochrony przed zanieczysz- 
czeniem wód morskich w pasie przybrzeżnym. Mając na względzie ogólno- 
krajową funkcję rekreacji i wypoczynku ludności, jaką spełnia pas nadmorski 
województwa szczecińskiego [2] przewiduje się budowę kilku średniej 
wielkości oczyszczalni ścieków gospodarczo-bytowych. Są to: biologiczne 
oczyszczalnie ścieków w Świnoujściu i Międzyzdrojach, mechaniczno-bio- 
logiczne oczyszczalnie dla Mrzeżyna i Trzebiatowa, grupowe oczyszczalnie 
mechaniczno-biologiczne w Kołczewie, dla Dziwnowa-Dziwnówka oraz 
zakończenie zbiorczej oczyszczalni ścieków w Pobierowie. Budowa wymie- 
nionych oczyszczalni w miejscowościach nadmorskich jest niezbędna dla 
utrzymania właściwych warunków sanitarnych ludności sezonowej, zakłóca- 
nych obecnie koniecznością wywozu ścieków i osadów transportem samo- 
chodowym. 

Przedstawione zamierzenia budowy komunalnych i przemysłowych oczysz- 
czalni ścieków uwzględniają zadania, wynikające z Konwencji o Ochronie 
Morskiego Środowiska Obszaru Morza Bałtyckiego z 1974 r. (Konwencja 
Helsińska), której sygnatariuszem jest Polska. Zadania te obejmowały między 
innymi na terenie województwa szczecińskiego budowy oczyszczalni ścieków 
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w Szczecinie, Świnoujściu, Międzyzdrojach, Wolinie, Pobierowie i w Dziw- 
nowie. 

Ważnym elementem ochrony wód przed zanieczyszczeniem, w sytuacji 
stałego zmniejszania się zasobów wody wysokiej jakości, jest wymienione 
w programie [1] realizowanie stref ochronnych w rejonach posiadających 
takie zasoby. Strefy ochronne zasobów wodnych, nadających się do celów 
wodociągowych będą realizowane w województwie szczecińskim w rejonie 
aglomeracji (Mścięcino, Międzyodrze) oraz w rejonie pasa nadmorskiego 
(zlewnie rzek Wołczenicy i Grzybnicy wraz z jeziorem Ostrowo). Będą 
również finałizowane prace w strefie ochronnej jeziora Miedwie, zasilającego 
wodociąg szczeciński. 

Program unieszkodliwiania i utylizacji odpadów bytowo-gospodarczych 
[1] przewiduje rozpoczęcie w 1995 r. budowy kompostowni lub spalarni 
odpadów w dziewięciu największych aglomeracjach kraju, w tym w aglomera- 
cji szczecińskiej. 

Ze względu na wyczerpanie w ciągu najbliższych lat możliwości składowa- 
nia odpadów gospodarczo-bytowych na trzech istniejących wysypiskach aglo- 
meracji szczecińskiej (przy ul. Podbórzańskiej, przy ul. Komety i w Dołujach) 
i wyłączenie ich z eksploatacji, odpady składowane będą na jednym składowis- 
ku w Sierakowie. Jego pojemność (po powiększeniu) wystarczy zaledwie do 
1995 r. Stąd występuje konieczność rozpoczęcia budowy kompostowni 
w Szczecinie o wydatku 360 t/r w 1989 r. tak, aby w 1995 r. rozpoczęła 
produkcję. 

Do odpadów szczególnie niebezpiecznych dla środowiska, poza wymienio- 
nymi w programie [1] (pestycydy, odpady porafineryjne, koksownicze i lakier- 
nicze) należy zaliczyć fosfogipsy. Ze względu na ich ilość (str. 6) oraz toksyczne 
oddziaływanie na wody powierzchniowe i podziemne, problemy technologii 
gospodarczego wykorzystania oraz bezpiecznego składowania fosfogipsów 
powinny być ujęte w programie narodowym. 

W programie ochrony powierzchni ziemi [1] nie uwzględniony został 
problem ochrony brzegów morskich przed procesami abrazji. Długość brzegu 
morskiego wraz z brzegami morskimi wód wewnętrznych wynosi w wojewódz- 
twie szczecińskim 323,5 km, w tym 241 km brzegów Zalewu Szczecińskiego. 
W wyniku stałego wzrostu procesów abrazyjnych, program ochrony brzegów 
morskich uwzględnia stały wzrost robót ochrony technicznej i biotechnicznej. 
Program ochrony brzegów Zalewu Szczecińskiego uwzględnia stałą konserwa- 
cję i odbudowę obwałowań przeciwpowodziowych (o długości 88 km) w wyni- 
ku procesów osiadania słabonośnego podłoża gruntowego i zanikania natural- 
nych porostów trzciny, hamujących działania abrazyjne. 

Nie przedstawiono w programie [1] podstawowych dla gospodarki wodą, 
wyników bilansowania potrzeb i zasobów wody dla zaopatrzenia gospodarki 
komunalnej (woda wysokiej jakości) przemysłu i rolnictwa w roku 2010. 
Przedstawiona informacja o sumarycznym wzroście potrzeb wody z 15,6 
mld m3 w roku 1986, o 9,5 mld m*/r według wariantu A (program potrzeb) 
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albo o 8,6 mld m*/r (program możliwości) w roku 2010 nie obejmuje potrzeb 
wodociągów wiejskich oraz rolnictwa. Zarejestrowanych poborów wody w ro- 
ku 1986 nie można utożsamiać z rozmiarami potrzeb, bowiem nie wszystkie 
potrzeby mogły być zaspokojone, wobec występowania znaczących deficytów 
wody. 

Synteza programu likwidacji deficytów wody w skali kraju nie została 
przedstawiona w programie [1] w sposób dostateczny. Przedstawione na 
mapie [1] proponowane przerzuty wody i kierunki zaopatrzenia w wodę nie są 
kompletne. Pominięto przewidywany przerzut wody o długości ponad 40 km 
ze zlewni rzek Wołczenicy i Grzybnicy (jezioro Ostrowo), w celu zaopatrzenia 
w wodę Świnoujścia i 14 miejscowości uzdrowiskowych w pasie nadmorskim 
województwa szczecińskiego. Konieczność przerzutu wody wynika z braku 
zasobów „słodkiej” wody na wyspach Uznam i Wolin oraz potrzeby zaopat- 
rzenia w wodę około 5,5 mln ludności z całego kraju, przebywającej w okresie 
sezonu w pasie nadmorskim województwa szczecińskiego. 

W programie narodowym pominięta została również rozpoczęta budowa 
drugiego rurociągu magistrali wodociągowej o długości 27 km, łączącej stację 
wodociągową nad jeziorem Miedwie ze Szczecinem. Realizacja drugiego 
rurociągu magistrali umożliwi zwiększenie wydajności istniejącego wodociągu 
do 200 tys. m*/d, po wykonaniu projektowanego przerzutu wody z rzeki 
Drawy do rzeki Płoni. 

W opisie mechanizmów ekonomicznych, stymulujących podejmowanie 
przedsięwzięć w zakresie ochrony środowiska i gospodarki wodnej nie wymie- 
niono powoływania i działalności spółek wodno-ściekowych. Praktyczna 
działalność spółek, powstałych w województwie szczecińskim („Świna”, „Trze- 
biatów-Mrzeżyno”, „Parkowa”, „Północ”) potwierdziła celowość powoływania 
spółek oraz dalszego zwiększania ich efektywności. Ten sposób podejmowania 
przedsięwzięć inwestycyjnych, ich realizacja i użytkowanie wykonanych obiek- 
tów, okazał się szczególnie efektywny w warunkach użytkowania obiektów 
rekreacyjnych w pasie nadmorskim oraz na terenie aglomeracji. 


PODSUMOWANIE 


Porównanie oceny istniejącego stanu i koncepcji ochrony środowiska 
przyrodniczego, przedstawionych w programie narodowym, do występujących 
w województwie szczecińskim prowadzi do następujących wniosków: 

Stan istniejący środowiska przyrodniczego w województwie szczecińskim 
jest gorszy, niż to wynika z oceny przedstawionej w programie narodowym, 
zwłaszcza w odniesieniu do zanieczyszczenia powietrza atmosferycznego i wód 
rzeki Odry. 

W ocenie stanu istniejącego, a zwłaszcza w koncepcji ochrony środowiska 
przyrodniczego i gospodarki wodnej, pominięto szereg elementów spe- 
cyficznych dla województwa szczecińskiego i jego roli w gospodarczym 
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układzie kraju, jak gospodarka morsko-przemysłowa i funkcja turystyczna 
pasa nadmorskiego. 

Przedstawiony projekt narodowego programu ochrony środowiska przyro- 
dniczego nie jest w pełni odpowiednią podstawą do opracowania programu dla 
województwa szczecińskiego. 

W związku ze specyfiką regionu, potrzeba opracowania wydzielonego 
kompleksowego programu ochrony środowiska naturalnego i gospodarki 
wodnej województwa szczecińskiego jest w pełni uzasadniona. 


JERZY BOCZAR 


PROBLEMS OF SZCZECIN DISTRICT IN ACCORDING 
WITH THE NATIONAL ENVIRONMENTIAL PROTECTION 
UNTIL 2010 


Summary 


Some analysis convergence of assessement of the state of natural environment and of the water 
supply of Szczecin District has been shown in the paper. This analysis has been presented in 
environmential protection with the reconnaissance which has been done at that district. 

The general conception of environmential protection and water conservation, which has been 
presented in the national program, also has been compared to the principal directions of 
functioning operations. 

These operations have been elaborated by the professional specialists of Szczecin District. 


EXH BOUAP 


IIPOEBJIEMATAKA INENAHCKOTO BOEBOĄCTBA 
COOTHOCHTEJIGHO C HAPOJNHOM IIPOTPAMMOM 
OXPAHBI IIPAPOJHOM CPEJIbl JIO 2010 F. 


Pe3roMe 


Pazpa6oTaH aHaJIH3 CXOXĄCHAA OLNEHKA COCTOAHHA IPAPO1HOŃ CPENBI H CHAaÓXKEHHA BONO 
B II[eqAACKOM BOEBOĄCTBE, Npe1CTaBJIeHHOŃ B HapONHOŃ IporpaMMe, C OÓCJIeNOBAHHEM, 
HpoBe1eHHHIM B BOEBOJĄCTBe. CpaBHeHa KOHHENNIAA OXPAHKI IPAPOĄAOCŃH CPEĄBI HU BOJNHOTO 
XOBAICTBA, HpeNIoxeHHaA B HapONHOH HporpaMMe, C TJIaBAHKIMA HanpaBJIeHAAMA NeńcTBHii 
B 3TOH OGJacTu, pazpaGoTaHHBIMH KOJLIEKTABAMA CIELAAJIACTOB IIIenAACKOTO BOEBOJĄCTBA. 
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WPŁYW NIEMONOTONICZNYCH ZMIAN 
PEWNEJ KLASY FUNKCJI 
LEPKOŚCI CIECZY NEWTONOWSKIEJ 
NA ROZKŁAD PRĘDKOŚCI CIECZY 
W GŁĘBI PÓŁPRZESTRZENI 
OGRANICZONEJ PORUSZAJĄCĄ SIĘ PŁYTĄ 
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ul. Wielkopolska 15, 70-451 Szczecin; 
Katedra Techniki Cieplnej Politechniki Szczecińskiej, 
al. Piastów 19, 70-310 Szczecin 


Streszczenie. W pracy ukazano analizę wpływu niemonotonicznych zmian 
wartości funkcji lepkości cieczy newtonowskiej w głębi półprzestrzeni na rozkład 
prędkości cieczy wywołany przez poruszającą się ze stałą prędkością nieograniczoną 
płytę na brzegu tej półprzestrzeni. Z rozważań wynika, że ekstrema lokalne dla klasy 
wymiernych funkcji lepkości nie wywołują niemonotonicznych zmian funkcji roz- 
kładu prędkości cieczy w głębi półprzestrzeni. Wykazano, że dla tej klasy funkcje 
rozkładu prędkości są Ściśle malejące. 

Słowa kluczowe: Przepływ cieczy w mieszalniku — niemonotonicznie zmien- 
na lepkość cieczy — analityczne rozwiązanie problemu. 


WSTĘP 


Niniejsza praca stanowi próbę teoretycznego opisu przepływu cieczy, który 
występuje w trakcie jej mieszania. W prostokątnym układzie współrzędnych 
płyta o nieograniczonej długości porusza się ruchem płaskim i jednostajnym 
w kierunku osi x z prędkością dodatnią U =const w sąsiedztwie newtonowskiej 
lepkiej cieczy (ryc. 1). Przyjmujemy, że w chwili czasu £t=0 ciecz nie porusza się, 
a jej ruch rozpoczyna się dla t>0. Ponadto zakładamy, że gęstość cieczy p jest 
w badanym przepływie stała. Tak więc ciecz jest nieściśliwa, a jej przepływ 
rozumiemy jako izotermiczny, ustalony. 

Problemem niniejszej pracy jest analityczne wyznaczanie rozkładu prędko- 
ści cieczy w półprzestrzeni, tzn. w kierunku osi y, dla cieczy newtonowskiej 
o niemonotonicznie zmiennej lepkości w głębi tej półprzestrzeni. Należy 
wyraźnie podkreślić, że lepkość cieczy nie zależy ani od ciśnienia, ani od 
temperatury, ponieważ przyjmujemy stałe wartości tych wielkości. Lepkość 
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ółprzestrzeni 


Ż.. 


ść pó 


głęboko 


nieograniczona długość płyty 


Rys. 1. Ruch nieograniczonej płyty o stałej prędkości U w sąsiedztwie półprzestrzeni wypełnionej 
newtonowską lepką cieczą o prędkości V,=V,(y, t); 1) lepka newtonowska ciecz, 2) linie prądu 
cieczy, 3) płyta 


zmienia się wyłącznie wraz ze zmianą składników chemicznych rozpatrywa- 
nych mieszanych cieczy. Ponadto lepkość zależeć będzie od czasu, co jest 
wywołane starzeniem się cieczy w trakcie jej przepływu. 

Nowym aspektem rozpatrywanego zagadnienia jest analiza wpływu niemo- 
notonicznej zmienności współczynnika lepkości kinematycznej w kierunku 
prostopadłym do poruszającej się płyty na wartości rozkładu składowej 
prędkości cieczy w tym samym kierunku. Szczególny przypadek przedstawio- 
nego problemu, ale dla cieczy o stałej lepkości, jest opracowany w literaturze 
mechaniki płynów [2]. 


WYPROWADZENIE RÓWNANIA 


Zwróćmy uwagę, że w rozpatrywanym problemie współrzędne V,, V,, V, 
wektora prędkości 6 cieczy są następujące: 
K=Kv.d. WV.d=0, K(y, =0, (1) 
gdzie O£y<+w, 0<xt<+w. 


Ponieważ przepływ cieczy wywołany jest ruchem płaskiej płyty, wobec tego 
naprężenia normalne w płynie nie występują, stąd 


p=const, 2) 


gdzie p oznacza ciśnienie. 
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Przyjmujemy, że lepkość kinematyczna v zmienia się po głębokości 
półprzestrzeni i czasie, zatem 


v=v(y, t). (3) 


Biorąc pod uwagę fakt, iż p=const oraz przyjęte założenia, równania ruchu 
upraszczają się do następującej postaci [2], (patrz Appendix 1): 


0%_0[,,,98% 
dr ay] 90 dy | (4) 


gdzie niewiadomą jest funkcja V,=V,(y, t). 
Zauważmy, że dla t=0 ciecz nie porusza się; natomiast jej ruch rozpoczyna 
się dla t>0, stąd mamy: 


V,v=0, t>0)=U, 6) 


gdzie U oznacza prędkość stałą. 
Z, drugiej strony w znacznej odległości od płyty, tj. dla y> + oo, ruch cieczy 
zanika. Zatem 


V,.W> +0, t>0)=0. (6) 
Zależności (5) oraz (6) są warunkami brzegowymi dla równania (4). 


CAŁKOWANIE RÓWNANIA RUCHU CIECZY NEWTONOWSKIEJ 
O ZMIENNEJ LEPKOŚCI 


W celu rozwiązania równania (4) zakładamy, że lepkość kinematyczna v jest 
rozpatrywana w następującej klasie funkcji: 


vy,d=v(ć), (7) 
gdzie bezwymiarowy argument 6 ma postać 
—— 8 
2/vt ©) 
Występujący we wzorze (8) charakterystyczny współczynnik wymiarowych 
wartości lepkości kinematycznej vy jest zdefiniowany wzorem: 


ć(, )= 


vy, t) 
Yyy = , 9 
, v+(v, t) ( ) 
gdzie v, jest funkcją określającą lepkość kinematyczną o wartościach bez- 


wymiarowych. 

Biorąc pod uwagę funkcję określoną wzorem (8) łatwo spostrzegamy, że 
gdyt>0iy>0,toć>+o;gdy0<t<+oiy—0,toć —0; natomiast gdy 
0<i<+o i y> +o, to 6 > +w. Ponadto o rozpatrywanej cieczy za- 
kładamy, że jej lepkość przed rozpoczęciem ruchu jest równa lepkości 
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w warstwie przyściennej po jego rozpoczęciu. Stąd na mocy uzyskanych wyżej 
granic oraz założeń (7), (8) i (9), mamy: 


M 
(FOJ=VO0=". (10) 
[) 
gdzie M jest dowolną stałą dodatnią. 
Aby rozwiązać równanie (4), przyjmujemy następujące podstawienie: 


V.(6)=U-/(ó), (11) 


gdzie zmienną € określa wzór (8), natomiast f(6ć) oznacza nieznaną funkcję 
o wartościach bezwymiarowych. 
Równanie (4) przedstawiamy teraz w następującej postaci: 
ÓV, óv0V,  GV, 
=— +v 2 * 
Ot dy dy dy 


Po wykonaniu różniczkowania równań (7), (8) oraz (11) otrzymujemy, [1]: 


(12) 


dV,_y 9%. 
W 227 (13) 
dv_dvóć 14 
dy dć dy” (14) 
0V, „,df 06 
0 dz dy (15) 
oraz 
8%, ..dYfo0y daf o 
aż" = (0) dĘ dy?” > 
gdzie 
ŻĘ. 1 Eo 6R: 1, 
ów zj O ŚR ŻE (17) 


Po wprowadzeniu zależności (13) — (17) do równania (12) oraz po uwzględ- 
nieniu wzoru (9) i po przeprowadzeniu odpowiednich przekształceń, otrzymuje- 


my: 
dźf_ 1 df/dy, 
wać (gę+2:). (18) 


W świetle podstawień (8) oraz (11) określamy teraz warunki brzegowe dla 
równania (18). Zauważmy więc, że jeśli y=0it>0, to 4=0. Wówczas na mocy 
(5) i (11), mamy: 


fE=0=1. (13) 
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Natomiast jeśli y= +00 i t>0, to Ę > + o. Wtedy na mocy (6) oraz (11), 
mamy: 
fE> +0)=0. (20) 
Równanie (18) po rozdzieleniu zmiennych i po jednokrotnym całkowaniu 
otrzymuje następującą postać: 


5 
daf__ Ci | 26: | 
„EEE OB na | GRE 
dę vlć) 7 | „(6 (21) 
[0 
gdzie v,(ć,)40 dla O<ć,<6<+w, a C; jest stałą całkowania. 


Po wykonaniu całkowania równania (21) otrzymujemy następujące roz- 
wiązanie ogólne równania (18), [3]: 


1 26 
J() -0 | exp |- Ę V "ĆJ a, dć, + C; , (22) 
[o 


gdzie v,(6,)+0 dla OKE, KE, <E<+w, a C, jest stałą całkowania. 
Zatem po nałożeniu na rozwiązanie ogólne (22) warunków brzegowych (19) 
i (20), uzyskujemy: 


+ 


oo Ś2 
I Ć sę 
=— R OZ d d =L1. 
Ci ||-ge|- I wić) | | AE (23) 
[0) 


Po uwzględnieniu w rozwiązaniu ogólnym (22) wzorów (23), a następnie po 
wykonaniu odpowiednich przekształceń oraz po wykorzystaniu zależności (8) 
i (11) uzyskujemy następującą ostateczną postać rozwiązania szczególnego 


równania (4): 
y 
V,(6)=V,|—— = 
AG Gz) 


y 52 j 
2yvot | PALĄ 
d d 
| V4 (Ć2) u |- [6 V+ (6,) a | ś2 


=U<1—9 


(24) 


+0 , 


1 26, 
| „() *? |- Ę: mJ de, de 


[U 


gdzie 0<6, <6, <6= wzwśkih 

vt 
W szczególnym przypadku, gdy funkcja lepkości kinematycznej v, o wartoś- 
ciach bezwymiarowych jest stała i równa jedności, to uwzględniając całkę 
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Eulera-Poissona, uzyskujemy [1]: 


y 2 (złe s 
"KIF |= z —ć)dć, | , 
wsk z)” | p| C9 J 05) 


gdzie 0£Ę,<$=" Naaźc, 


2./v5t 
Rozwiązanie (25) jest znanym rozwiązaniem Rayleigh'a w omawianym 
problemie dla cieczy newtonowskiej, ale o stałej lepkości kinematycznej. 


ANALIZA ROZWIĄZAŃ RÓWNANIA 
RUCHU CIECZY NEWTONOWSKIEJ DLA PEWNEJ KLASY 
NIEMONOTONICZNYCH FUNKCJI LEPKOŚCI 


Obieramy funkcję v, =v, (6) tak, aby bez utraty ogólności dalszych roz- 
ważań, była spełniona zależność (10), w której M/v,=1. Przyjmujemy zatem, że 


=zzż=ct |, (26) 


gdzie O£ć<+o, O0<b<+0 oraz —2/b<a<+o. 

Po zbadaniu zmienności funkcji (26) współczynnik a został tak dobrany, 

aby v, (6)>0. Kilka wykresów funkcji v, =v, (6) z badanej klasy (26) ukazano 
na ryc. 2. 
Funkcje v,=v,(ć) posiadają ekstremum lokalne dla ć=,/b, a mianowicie: 
maksimum lokalne dla 0 <a+ oo , minimum lokalne dla —2,/b<a<0. W tych 
przypadkach istnieje dla funkcji klasy (26) asymptota pozioma o równaniu 
v,=1. Natomiast gdy. a=0, to v,(ć6)=1 i wówczas lepkość jest stała. 
W wypadku gdyby a= —2/b, to odpowiednia krzywa na ryc. 2 oznaczałaby 
taki niemonotoniczny rozkład lepkości kinematycznej, dla którego w punkcie 
$=„/b lepkość dynamiczna cieczy „ dążyłaby do zera, a gęstość p miałaby 
wartość skończoną, ponieważ v=u/p. 

Po podstawieniu wzoru klasy funkcji (26) do rozwiązania szczególnego (24) 
równania (4) oraz po wykonaniu odpowiednich przekształceń, otrzymujemy: 


Ś 2 
J $(ć.) exp [-[ |śi $(ć,) dć,] dć, 
V,()=U41-0— CH ! (27) 
J (3) exp [— [ót 9(E,) dć,] dć, 
0 0 


gdzie 
6 +b 


pa ŚR 2, idą, dE5 
a gt i, 2 (28) 


$(ć)= 
y 


oraz 0£ĘKE=——=< tw. 
yt 
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bezwymiarowa lepkość cieczy 


0 1 l, 7 EE. 
bezwymiarowa głębokość półprzestrzeni  |P ŻWtb 


ue 
Ryc. 2. Rozkład bezwymiarowych wartości v, FE wymiernych funkcji lepkości kinematycz- 
+ 


nej w głąb półprzestrzeni dla O<ę*=ć//b<+oo oraz O<b<+w 
Rozpatrujemy teraz znak wyróżnika 4=a*—4b trójmianu występującego 
we wzorze (28). Zatem wzór funkcji (27) przybiera odpowiednie postaci dla 
4 <0, dla 4=0 i dla 4>0. Stąd funkcję zdefiniowaną wzorem (27) oznaczamy 
odpowiednio symbolami: V.qy, Vxqay 1 Va): W poniżej rozpatrywanych 
przypadkach uwzględniamy fakt, iż —2/b<a<+w. 


Przypadek 1 
Jeżeli 4 <0, to —2/b<a<24/b i0<b<+oo. Wtedy istnieje w punkcie 


ć=vb dla klasy funkcji (26) minimum lokalne gdy —2/b<a<0, natomiast 


maksimum lokalne gdy 0 <a<2/b. Po wyznaczeniu odpowiednich całek we 
wzorze (27), otrzymujemy: 


$ 
J A(€2) exp [B(6)+ C(62)] dź, 
Vw6©=Uil-z 


| A(2) exp [B(E,)+C(€,)] dć, 
gdzie 


6ż+b 


08 Prag 
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Ly, 2 —2b SNIDEZAW 
PROZĘ arctg w=uł (31) 


C(Ę,)= —63+2ać,, (32) 


dla —2/b<a<2/b, 0<b<+w i O<t,źt<+o. 


Zwróćmy uwagę na to, że gdy a—0, to A(6,)1, B(Ę,)>1 
i C(ę,)=—ćż. Wtedy otrzymujemy szczególny przypadek opisany wzorem 
(25). 
Przypadek 2 


Jeżeli 4=0, to a=2/b i0<b< + oo. Wówczas istnieje dla klasy funkcji 


(26) maksimum lokalne w punkcie c=4/b. Wtedy po wyznaczeniu odpowied- 
nich całek we wzorze (27), mamy: 


Ś 
J D(62) exp [E(62)] deą 
V2) (6)=U41— RZ ? (33) 
| D(6,) exp [E(6,)] dó, 
gdzie 
62 +b 
D PE SEO (34) 
(ć2) TEWSCE 
4b/b 


E(6,)=—63+4./be, — 


GŁ b Sa 


dla a=2/b, 0<b<+0 i O<Ę,<Ę<+w. 
Ryc. 3. przedstawia przebieg zmienności rozkładu prędkości V,,=V.2, (6). 
Przypadek 3 


Jeśli 4>0, to 2/b<a<+0 i O0<b<+o. Wtedy istnieje dla klasy 


funkcji (26) maksimum lokalne w punkcie ć =yb. Zatem po wyznaczeniu 
odpowiednich całek we wzorze (27) otrzymuje on następującą postać: 


$ 
[LF (ć2)”" $(6,) exp [C(62)] dć, 

V, a, (6)= U 1-55 > (36) 
J LF(E1?" $ (62) exp [C(62)] dć> 


gdzie 
F(ć)=(2—ćm) *' (62—60**. (37) 
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= 
E 


c| 


a=slb dla s=2 i b=1;2. 


| 
Va) 


U 


bezwymiarowa prędkość cieczy 


© 


10 20 AT 
bezwymiarowa głębokość półprzestrzeni 


Ryc. 3. Porównanie rozkładu prędkości cieczy newtonowskiej w głębi półprzestrzeni; 1) rozkład 

prędkości Rayleigh'a, 2) rozkład prędkości cieczy dla niemonotonicznej funkcji lepkości gdy b=1, 

3) rozkład prędkości cieczy dla niemonotonicznej funkcji lepkości gdy b =2; rozkłady 2 i 3 opisane 
są funkcją określoną wzorem (33) 


oraz 
ri=fifg dh 
paz GUS ERZE dla i=1, 2; (38) 
2 /a?—4b 
ć1=2(—a—4/a*—4b)<0, 6n=3(-a+./a?—4b)<0; (39) 


dla 2/b<a< +0, 0<b<+wo i OSĘ,<ć6<+w0; natomiast funkcje 
$=$(ć,) i C=C(6,) określają wzory (28) i (32). 

Po zbadaniu własności rozwiązania szczególnego (24) równania (4), w któ- 
rym funkcja lepkości kinematycznej jest określona wzorem (26) możemy 
sformułować następujące 


Twierdzenie 1 

Jeżeli istnieje klasa funkcji lepkości kinematycznej cieczy newtonowskiej 
o identycznych dodatnich wartościach w nieskończoności, równym wartościom 
na brzegu półprzestrzeni i posiada w jej głębi dodatnie ekstremum lokalne, to 
rozkład składowej prędkości cieczy wywołany przez poruszającą się ze stałą 
prędkością płytę na brzegu tej półprzestrzeni jest w jej głębi funkcją ściśle 
malejącą. 


4 Szczecińskie Rocz. Nauk. t. VI, z 1 
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Dowód 

1. Niech klasa funkcji lepkości będzie przyjętą klasą (26). 

Zauważmy zatem, że gdy 6 > +o, to v,(6)—1 i również gdy Ę—0, 
to v,(ć)—1. Ponadto dla O<ć<+w i —2/b<a<+0o pochodne klasy 
funkcji (26) określone wzorem: 


dE 7 —a) (6? -b)(6*+b)"*, (40) 
posiadają w badanym przedziale miejsce zerowe dla c=4b. Ponadto funkcje 


klasy (40) zmieniają znak w otoczeniu punktu ć =4/b oraz przy przyjętych 
założeniach 


D<<=+1>6; 
v. (Yb) NIM (41) 
Oznacza to, iż funkcje klasy (26) posiadają dodatnie ekstremum w punkcie 
ć=y/b. A zatem istnieje klasa funkcji lepkości kinematycznej spełniająca 
założenia twierdzenia 1. 
2. W celu wykazania ścisłej monotoniczności funkcji V,,,„,=V.,,(ć) dla 
j=1,2,3 określonych wzorami (29), (33) i (36) zanalizujemy teraz pochodne 
tych funkcji, które mają następujące postaci: 


Myy 460 «wp [B(0+CE] = 
d + o , 
; | AG) ep [BE)+CE] dt, 
Wa__y__ DOepIE(G] | i 
z | De)esp CE(EJJdE, 
dv, [F(G?6(E Cc 
U LF(6)* 9 (6) exp [C(6)] : (44) 


I [FEI 6(6) exp [C(EJ] de, 


[U 


gdzie U jest stałą dodatnią. 

Zauważmy najpierw, że funkcje A=A(6), D=D(ć), F=F(ć) i $=4(6) 
określone wzorami (30), (34), (36) i (28) dla O<ć,<Ę<+ oo nie posiadają 
miejsc zerowych. Istniejące zera funkcji (30), (34) i (36) nie należą do badanego 
przedziału. Łatwo spostrzec, że funkcje A=4 (6), D=D(), F=F(ć6)i$=4 (6) 
występują we wzorach (42),(43) i (44) w postaci czynników, przy czym 
pozostałe funkcje nie zerują się. Stąd pochodne podane wzorami (42), (43) i (44) 
również nie mają miejsc zerowych. Oznacza to, że badane funkcje (29), (33) 
i (36) nie posiadają ekstremum w badanym przedziale. 
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Zanalizujemy teraz znak poszczególnych czynników występujących we 
wzorach funkcji (42), (43) i (44). Wykorzystując podstawowe własności funkcji 
kwadratowej łatwo stwierdzamy, że funkcje A=A(6), D=D(ć) i F=F(ć) mają 
w badanym przedziale dodatnie wartości. Zatem wszystkie czynniki we 
wzorach (42), (43) i (44) są dodatnie z wyjątkiem czynnika (— U), bo czynnik 
U jest dodatni. Wobec tego badane pochodne (42) — (44) mają dla 
0x6, <6< + © wartości ujemne i nie zerują się. Oznacza to, że dla przyjętej 
klasy funkcji (26) funkcje V.vy=Kut), Va=Tx2(ó) I Ka "Va (6) 
określone wzorami (29), (30) i (36) są funkcjami ściśle malejącymi, co kończy 
dowód twierdzenia 1. 


ZASTOSOWANIE 


Rozpatrywany w niniejszej pracy problem znajduje zastosowanie m.in. 
w mieszalnictwie, a szczególnie tam gdzie stosuje się mieszadła turbinowe 
i łapowe [4]. W zbiorniku ulegają mieszaniu rozmaite ciecze o różnych 
lepkościach. W związku z tym ich mieszanina charakteryzuje się niemonotoni- 
cznie zmienną lepkością w przestrzeni między poruszającym się mieszadłem 
a głębią cieczy. Zatem problem rozkładu prędkości cieczy w rozpatrywanej 
półprzestrzeni ma istotne znaczenie w zagadnieniach mieszalnictwa (ryc. 4). 


Ryc. 4. Mieszadło z pionowymi płytkami o ruchu posuwisto-zwrotnym. Zakłada się, że: U;*0, 
T=const, p=const, V,=V,(y, t), V,=0, V,=0 


4* 
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WNIOSKI 


Zasadniczym wynikiem rozważań przedstawionych w niniejszej pracy jest 
analiza wpływu niemonotonicznych zmian funkcji lepkości kinematycznej 
cieczy newtonowskiej w głębi półprzestrzeni na rozkład prędkości cieczy 
wywołany przez poruszającą się ze stałą prędkością nieograniczoną płytę na 
brzegu tej półprzestrzeni. 

Z przedstawionych badań wynika, że ekstrema lokalne dla klasy wymier- 
nych funkcji lepkości kinematycznej nie wywołują niemonotonicznych zmian 
funkcji rozkładu składowej prędkości cieczy w głębi półprzestrzeni. 

Uzyskane w niniejszej pracy funkcje rozkładu prędkości cieczy, dla przyję- 
tej klasy wymiernych funkcji lepkości, są ściśle malejące — podobnie jak 
w klasycznym problemie Rayleigh'a. Fakt ten nie oznacza jednak, że mono- 
toniczny spadek wartości funkcji składowej prędkości cieczy newtonowskiej 
o niemonotonicznym zmiennym współczynniku lepkości jest identyczny z mo- 
notonicznym spadkiem wartości funkcji składowej prędkości cieczy w głębi 
półprzestrzeni dla klasycznego problemu Rayleigha. 


Appendix 1 


Problem pracy opisany jest równaniami pędu, które dla stacjonarnego 
laminarnego przepływu mają postać [2]: 


dv, __0p 8 a 
Pac gz tz CAO p 
dV, 6p 3 
Zakł ZEREM, Ka 
Pt GE (2u6,,)+  210,)+ (2 1O,.), (A.1) 
SĄ 210,,)+>- (210,,)+> (2u6,, 
P dt 82” dx 218.,)+ z, | » Duryk LQ.) 


oraz równania ciągłości o postaci [2] 


ÓV, ÓW, OV, 


-0 
0 0 2 3 


(A.2) 


gdzie współrzędne tensora prędkości deformacji są następujące [2]: 


0V, dY, 0V, 
0. 0x” „= Za 
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_ o _1/(0V, OV, (A.3) 
B= B= (7 * sk 
1/6V, dV, 
8,-0,—; (gat 37). 


a symbol u oznacza współczynnik lepkości dynamicznej. 
Wykorzystując założenia (1) i (2) oraz zależności (A.3) w równaniach (A.1) 
oraz (A.2) uzyskujemy postać równania (4). 


w 


SAWĘONOA 
1 op oo 
= kę 


RY 


x(1): 2)» V2) 


AGR ru m mm m O 
= 


ZESTAWIENIE OZNACZEŃ 


funkcje zmiennej ć, o wartościach bezwymiarowych 
funkcja zmiennej ć, o wartościach bezwymiarowych 

stałe całkowania o wartościach bezwymiarowych 

funkcje zmiennej ć, o wartościach bezwymiarowych 
dodatnia stała wymiarowa [m?/s] 

stałe bezwymiarowe 

prędkość płyty [m/s] 

współrzędne wektora prędkości cieczy po kierunkach x, y, z 
[m/s] | 

współrzędne wektora prędkości cieczy po kierunku x [m/s] 
liczba rzeczywista bezwymiarowa z przedziału (— 2,/b,+ co) 
liczba rzeczywista bezwymiarowa z przedziału (0,+ oo) 
funkcja zmiennej 6 o wartościach bczwymiarowych 
indeks=l1, 2 

indeks=1, 2, 3 

stała wartość ciśnienia [ Pa] 

czas [s] 

współrzędne układu kartezjańskiego [m] 

vektor prędkości cieczy [m/s] 

wyróżnik trójmianu kwadratowego o wartości bezwymiaro- 
wej 

współrzędne tensora prędkości deformacji [1/s] 

funkcja zmiennej ć, o wartościach bezwymiarowych 
współczynnik lepkości dynamicznej [Pa:s] 

lepkość kinematyczna płynu [m?/s] 

charakterystyczny współczynnik lepkości kinematycznej pły- 
nu [m?/s] 

lepkość kinematyczna płynu o wartościach bezwymiarowych 
bezwymiarowe argumenty (ć*=ć/,/b) 

bezwymiarowe argumenty jako parametry funkcji podcał- 
kowych 

bezwymiarowe ujemne miejsca zerowe trójmianu kwadrato- 
wego 

gęstość cieczy [kg/m*] 
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ANDRZEJ ANTONI CZAJKOWSKI, KRZYSZTOF WIERZCHOLSKI 


THE INFLUENCE OF NON-MONOTONE VALUES 
OF THE SOME CLASS OF VISCOSITY FUNCTION 
ON THE NEWTONIAN FLUID VELOCITY DISTRIBUTION 
IN THE HALF-SPACE RESTRICTED 
THROUGH THE MOVABLE PLATE 


Summary 


Some analysis of the influence of non-monotone values of the some class of viscosity function 
on the Newtonian fluid velocity distribution in the half-space restricted through the movable 
infinite plate with constant velocity has been given in the paper. 

Has been proved that the local extremum of the some class of the rational viscosity function of 
Newtonian fluid do not produce the non-monotone values of fluid velocity in the half-space. The 
fluid velocity values are strictly decreasing in the half-space deep. 


AH/DK3M AHTOHH UAAKOBCKA, KILŁIIITOG BEXXXOJIECKA 


BJIAAHAE HEMOHOTOHHBIX HU3MEHEHHIAi 
3HAUEHAA DYHKUMA BA3KOCTH 
HBIOTOHOBCKOŃ XHNKOCTH OIIPENEJJEHHOTO KJIACCA 
HA PACIIPEJIEJJEHAE CKOPOCTH XHHKOCTH 
B TJIYBHHE IIOJIYHIPOCTPAHCTBA, 
OFPAHHMUEHHOTO J|IBIOKYIIEACA IJIATOA 


Pe3sirome 


B paóore npeącTaBJieH AHAJIA3 BIAAHHA HeMOHOTOHHbIX H3MeHeHHi 3HaueHui byHKuAń 
BA3KOCTH HbIOTOHOBCKOŃ XHAJIKOCTH B TJIyÓHHE IOJIYIPOCTPARCTBAa Ha pacnpezreJleHAe CKOPOCTA 
XKHJĄKOCTH, BBI3BAHHOe€ HeOTpAHHUEHHOH IIHTOŃ, ABAKYlNeiCA C NOCTOAHHOI CKOPOCTBIO, Ha 
Kpało 3TOFO IOJIylMpocTpaHCcTBa. 

M3 paccyxNeHHf CJIENYET, UTO JIOKAJIBHBIE 3KCTDEMYMHI JU1A KJIACCA pALAOHAJIEHKMX PyYHKLUAŃ 
BA3KOCTH He BbI3BIBAIOT HEMOHOTOHHBIX H3MCHeHHii QyHKLUAŃi pacHpeqrejeHAA CKOpOCTH 
XKHJĄKOCTA B TJIYÓHHE NOJLyHPOCTPAHCTBA. 

YCTAaHOBJIEHO, 4TO MIA AAHHOTO KJIACca PyHKHHAH pacnpeNeJleAHA CKOPOCTA ABJIAŁOTCA ABHO 
yMEHBIIIAIOLNAMHCA. 
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JAN DUDCZAK 


SYNTEZA SYSTEMÓW PRZEMYSŁU CHEMICZNEGO 
ZE ZNANĄ RELACJĄ PORZĄDKUJĄCĄ PODSYSTEMY 


Instytut Inżynierii Chemicznej i Chemii Fizycznej Politechniki Szczecińskiej, 
al. Piastów 42, 71-065 Szczecin 


Streszczenie. Opracowano algorytm syntezy systemów przemysłu chemicz- 
nego, których podsystemy mogą zostać uporządkowane za pomocą pewnej relacji R. 
Algorytm oparty jest na strategii poszukiwawczej w dyskretnej przestrzeni stanów. 
Może on być zastosowany do syntezy systemów zasilanych pojedynczą surówką 
i wieloma strumieniami surówkowymi. Przedstawiono przykłady zastosowań al- 
gorytmu do syntezy instalacji oczyszczania gazów pochodzących z emisji wieloźród- 
łowej i syntezy węzła krystalizacji w instalacji regeneracji kąpieli koagulacyjnej 
w wytwórni włókien wiskozowych. 


Słowa kluczowe: synteza systemów — oczyszczanie gazów — krystalizacja. 


WSTĘP 


Synteza nowego systemu procesowego, w którym zadane strumienie 
surówkowe mają zostać przetworzone w żądane produkty w optymalny spo- 
sób, polega na określeniu: 

— optymalnej struktury systemu, tj. sposobu w jaki podsystemy powiązane są 
ze sobą i z otoczeniem systemu strumieniami procesowymi, 

— optymalnych wartości parametrów procesowych opisujących cechy podsys- 
temów i strumieni. 

Problem syntezy może być sformułowany jako następujący problem op- 
tymalizacji dyskretno-ciągłej. 

Poszukuje się 


MIN f[V,H] (1) 
V,H 


w obecności ograniczeń: 
— modeli matematycznych podsystemów 


M*ÓT[VÓ]=0, jeIS, (2) 
— równań powiązań 


S,„„=RC[S,,, H], (3) 


wy? 


56 


— ograniczeń nierównościowych na zmienne procesowe 


Z[V]<0 (4) 
— warunków zupełności struktury 


m 


He. 6) 
W powyższych formułach (ryc. 1): 


Wa (G) W - (i) 
wa Swe, Sw SZW 


Ryc. 1. Schemat strumieni zewnętrznych podsystemu 


— zbiór zmiennych procesowych 


V=S„US„UC, (6) 
— zbiór cech opisujących strumienie wejściowe 
S.= U Sie, (7) 
je(ISUJ0) 


zbiór cech opisujących strumienie wyjściowe 


Sy = U SO . (8) 
je(ISUJ0) 
— zbiór cech podsystemów ! 
C= U CU : (9) 
je 


— zbiór struktur zupełnych, tj. struktur zapewniających spełnienie ograni- 
czeń (4) 


A=(A.), zesz. (10) 


Wielowymiarowość obu przestrzeni decyzyjnych, nieliniowość problemu 
i ogromna różnorodność systemów procesowych sprawiają, że opracowanie 
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jednego uniwersalnego algorytmu, efektywnego w stosunku do wszystkich 
typów systemów procesowych jest praktycznie niemożliwe, wobec czego 
algorytmy syntezy opracowywane są oddzielnie dla systemów poszczególnych 
typów [1L, 2]. 


STRATEGIA POSZUKIWAWCZA W PRZESTRZENI STANÓW 


Ogólna strategia rozwiązywania problemu syntezy (1) — (5), która jest 
podstawą algorytmu syntezy przedstawionego w niniejszej pracy, oparta jest na 
strategii poszukiwawczej w dyskretnej przestrzeni stanów. Strategia ta wywo- 
dzi się z teorii rozwiązywania zadań, gałęzi teorii sztucznej inteligencji 
zajmującej się rozwiązywaniem problemów dyskretnych [3]. 

Proces rozwiązywania rozpatrywany jest tam jako poruszanie się w dys- 
kretnej przestrzeni stanów zadania. Poszukiwanie zaczyna się w pewnym stanie 
startowym (początkowym) i kończy w momencie osiągnięcia pewnego stanu 
docelowego, będącego szukanym rozwiązaniem. Elementarny ruch od stanu do 
stanu wynika z działania na stan macierzysty pewnego operatora zwanego 
operatorem następstwa. Zbiór operatorów następstwa działających na dany 
stan nazywany jest operatorem gałęzienia. Stan początkowy i operator 
gałęzienia opisują niejawnie całą przestrzeń stanów. W danym momencie 
poszukiwań do gałęzienia kandyduje zazwyczaj kilka stanów. Wybór stanu do 
gałęzienia zależy od użytej techniki przeglądu. Może ona być ślepa: wszerz 
(First-In-First-Out) lub w głąb (Last-In-First-Out) — gdy o gałęzieniu decyduje 
jedynie wiek stanu, lub ukierunkowana — gdy o gałęzieniu stanu decyduje 
wartość pewnej funkcji sterującej związanej z każdym stanem. Aby zmniejszyć 
liczbę przeglądanych stanów stosuje się często tzw. kryteria odsiewające 
(odcinające). Kryteria te, zwykle heurystyki ściśle związane z konkretnym 
uwarunkowaniem problemu, zawężają przestrzeń poszukiwań, lecz zazwyczaj 
za cenę rezygnacji z gwarancji osiągnięcia rzeczywistego stanu docelowego. 

Analiza algorytmów syntezy przedstawionych w literaturze, jak również 
algorytmów z dziedziny programowania całkowitoliczbowego i metod siecio- 
wych wykazuje, że wiele z tych algorytmów wykorzystuje strategię po- 
szukiwawczą w przestrzeni stanów, chociaż w większości przypadków fakt ten 
nie jest wyraźnie zaznaczony przez autorów. 

Aby sformułować algorytm poszukiwawczy w dyskretnej przestrzeni sta- 
nów, należy: zdefiniować pojęcie stanu zadania, sformułować operator gałęzie- 
nia, określić stan początkowy, opisać cechy stanu docelowego, wybrać technikę 
poszukiwań, (opcjonalnie) sformułować kryteria odsiewające. 

Kluczową decyzją, rzutującą na wszystkie pozostałe, jest wybór definicji 
stanu zadania. Analiza algorytmów syntezy przedstawionych w literaturze 
doprowadziła do wniosku, że istnieją trzy ogólne sposoby definicji stanu: 
— struktura cząstkowa, tj. struktura, w której spełnienie ograniczeń (4) jest 
niemożliwe, 

H!-HS'+0, dla każdego zeSZ, (11) 
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— struktura zupełna 


HSeHA, 
— zbiór struktur zupełnych 


HS=(A. :zeZZcSZ). (12) 
Szczegółowa postać stanu musi zostać wybrana w najwłaściwszy sposób dla 
konkretnego, aktualnie rozwiązywanego, problemu. Wybór ten powinien 
zapewniać możliwość sformułowania efektywnego operatora gałęzienia, który 
wygeneruje przestrzeń stanów o umiarkowanej liczbie stanów, wśród których 
znajduje się stan docelowy. 
Ogólną postacią stanu stosowaną w niniejszej pracy jest struktura zupełna. 
Z każdym stanem związana jest wartość funkcji oceniającej 
f* [HS] =MIN f[V, H]ys. 
v (13) 
Stanem docelowym jest więc stan charakteryzujący się najmniejszą wartością 
funkcji oceniającej. 
Przy formułowaniu szczegółowej postaci stanu przyjęto następujące założe- 
nia: 
— znany jest zbiór IT typów podsystemów, z których syntezowany jest system, 
— każdy stan zawiera co najmniej po jednym podsystemie każdego typu, 
— znane są maksymalne ilości podsystemów poszczególnych typów, które mogą 
wejść w skład syntezowanego systemu. 
Definicja operatora gałęzienia oparta jest na założeniu, że znana jest pewna 
relacja R porządkująca etykiety typów podsystemów 


..mRn... 


Relacja ta może być określona w różny sposób: heurystycznie, za pomocą 
wzoru, itp. 


SYSTEMY ZE ZNANĄ RELACJĄ PORZĄDKUJĄCĄ 


Liczba strumieni w systemie zależy od liczby podsystemów i jest od niej 
wielokrotnie większa. Sprawia to, że liczba możliwych struktur rośnie bardzo 
gwałtownie z liczbą podsystemów. Dla uproszczenia problemu podzielmy 
strumienie materiałowe występujące w systemie na strumienie kluczowe 
i strumienie pomocnicze. Podział ten jest dość subiektywny i trudno jest 
sformułować precyzyjną ogólną regułę podziału dającą jednoznaczną od- 
powiedź, szczególnie dla bardziej skomplikowanych systemów. Ogólnie można 
stwierdzić, że do zbioru strumieni kluczowych będą zaliczane strumienie 
pochodne od zadanych strumieni surówkowych, przetwarzane w systemie 
w celu osiągnięcia żądanych produktów. Wszystkie pozostałe strumienie 
zakwalifikowane będą do zbioru strumieni pomocniczych. 

Zbiór S cech opisujących strumienie materiałowe systemu dzieli się więc na 
dwa rozłączne podzbiory: 
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G — zbiór cech opisujących strumienie kluczowe, 

U — zbiór cech strumieni pomocniczych. 

Także zbiór strukturalny H dzieli się na dwa rozłączne podzbiory: 

HG — opisujący powiązania podsystemów strumieniami kluczowymi, 

HU — opisujący powiązania strumieniami pomocniczymi. 

Przyjmijmy założenie, że powiązania podsystemów strumieniami pomocniczy- 
mi wynikają jednoznacznie ze sposobu powiązań podsystemów strumieniami 
kluczowymi. Zakłada się więc, że znana jest funkcja określająca wartości 
elementów zbioru HU na podstawie znanych wartości elementów zbioru HG. 
W dalszej części pracy problem syntezy dotyczyć więc będzie wyłącznie zbioru 
HG. Każdy strumień kluczowy systemu można zapisać jako <z,i>, gdzie 
z przedstawia etykietę źródła strumienia, zaś i etykietę podsystemu, do którego 
strumień jest kierowany. Etykieta źródła ma postać z=(j, st), gdzie j jest 
etykietą podsystemu, zaś st jest etykietą typu strumienia. Etykieta podsystemu 
zawiera w sobie etykietę typu podsystemu. Przyjmuje się, że w systemie mogą 
występować strumienie dwóch typów: r i s (mogą to być na przykład dwa różne 
stany skupienia). Biorąc pod uwagę relację porządkującą R, strumienie łączące 
podsystemy <(j, st),i> można podzielić na: 

— strumienie sekwencyjne, dla których m,Rm,, 

— strumienie zwrotne (recyrkulacyjne), dla których m,RRm,. 

Przyjmijmy, że systemy rozpatrywane w niniejszej pracy będą zbudowane 
z podsystemów o ogólnej strukturze łączących strumieni zewnętrznych przed- 
stawionej na ryc. 2. 


REG 


U z] 
ELA 


Ryc. 2. Schemat kluczowych strumieni łączących podsystemu 


Powyższe ustalenia pozwalają teraz na zdefiniowanie przestrzeni stanów 
zadania syntezy systemów zbudowanych z podsystemów uporządkowanych 
relacją R. 

Każdemu strumieniowi surówkowemu przypisano tzw. linię procesową 
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imitującą oddzielny ciąg procesowy, w którym strumień ten może być 
przetwarzany całkowicie lub częściowo. Etykieta podsystemu ma postać (I, m), 
gdzie I=1,2,..., Nf jest etykietą linii procesowej, w której znajduje się 
podsystem, zaś m € IT jest etykietą typu podsystemu. Etykieta źródła strumie- 
nia ma postać z=(l, m, st), a etykieta strumienia <(/, m, st), (k,n)>. 

Stan zadania został zdefiniowany jako system, w którym każdy strumień 
przerabiany jest kolejno w podsystemach wszystkich typów. 

Stanem startowym jest system zbudowany z N/f niezależnych linii proceso- 
wych, z których każda zawiera wszystkie podsystemy ze zbioru IT powiązane 
jedynie strumieniami sekwencyjnymi 


<(l, m, r), (l, n): mRn>. 


Istnieją dwa operatory gałęzienia: 1) operator S zmieniający sposób powiąza- 
nia podsystemów strumieniami sekwencyjnymi, 2) operator R zmieniający 
sposób powiązania podsystemów strumieniami recyrkulacyjnymi. Nazwijmy 
strumieniem poprzecznym strumień sekwencyjny 


<(l, m, r), (k,n):k=1,2,...,l-1; mRn>. 


Wprowadzenie do struktury strumienia poprzecznego powoduje, że w linii 
procesowej I przestają istnieć podsystemy o typach n takich, że mRRn. 
Oznaczmy zbiór etykiet typów podsystemów pozostałych w linii ! jako IRT,. 
Operator S działa jedynie na stany utworzone działaniem operatora S. Tworzy 
on z danego stanu o strumieniach poprzecznych 


AUO LITECYZACH , <(yp> Myip> r), (KNip> Ny)> 
jego następniki przez dodanie strumienia poprzecznego 
<(ntp+ 1> MNtp+1> 7)» (kNip+ 1» ANip+1)> 
spełniającego (ryc. 3): 
lnp+15lnpt1: ---, Nf, 
kuip+171; ...y Lyip> 
NyNip+ 1 € IRĘ 1> 
Myip+1 Rhyg+1- (ryc. 3). 
Operator R tworzy z danego stanu o strumieniach poprzecznych 
<(l, , My > r), (ki > n;)>, ... „<(vp> MNip > 1 (KNip > Nyjp)> 
i strumieniach recyrkulacyjnych 
<(lh , my , sty) , (ką , n1)>, +... <(lv; MN, Sty) , (KNv; Ny) > 
jego następniki przez dodanie strumienia recyrkulacyjnego 
<(lvb+ 15 Myb+ 1> Styg+ 1) (KNb+ 1» MNb+ 1)> 
spełniającego 
INb+1=lnb> +» Nf, 
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IT=(A,B,C) 
ARBRCc 


—4A)3.8)A3.0>— Nf=3 


Ryc. 3. Przykład działania operatora gałęzienia S 
kyr+1=1, >, e» Nf, 
Myb+1 € (IRTJ,,,, i (my, RR myb+1, gdy lws+1 Int)» 


, LR 
Nyy+ 1 RRMyb+1- 


ZASTOSOWANIA ALGORYTMU 


Opracowany algorytm poszukiwawczy zastosowano do syntezy kilku 
wybranych systemów procesowych. Poniżej przedstawione będą wyniki dwu 
takich prac. Dokładny opis problemów można znaleźć w cytowanej literaturze. 


SYNTEZA INSTALACJI OCZYSZCZANIA GAZÓW [4] 


Do instalacji wchodzą trzy strumienie gazowe zanieczyszczone parami 
rozpuszczalnika organicznego. Cechy strumieni są następujące: 

F,=0,2 kmol/s, YF1 0,08, 

F,=0,075 kmol/s, YF2 = 0,02, 

F.=0,025 kmol/s, Ya =0,004. 

Temperatury i ciśnienia wszystkich strumieni są w przybliżeniu jednakowe 
i wynoszą odpowiednio 20”C i 10”N/m?. Ze strumieni tych należy usunąć 
składnik zanieczyszczający, tak, by jego stężenie w strumieniach opuszczają- 
cych instalację nie przekraczało 0,01% . Należy zsyntezować instalację o naj- 
niższym koszcie oczyszczania. 

Do oczyszczania gazu mogą być użyte trzy typy podsystemów oczysz- 
czających: niskotemperaturowa kondensacja (K), absorpcja (AB), adsorpcja 
(AD). Na podstawie danych literaturowych uporządkowano etykiety typów 
podsystemów w następujący sposób: 

K RABRAD. 
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Schemat łączących strumieni zewnętrznych podsystemu oczyszczającego cha- 
rakteryzuje się, w porównaniu z przyjętym uprzednio schematem ogólnym 
(ryc. 2), brakiem strumieni recyrkulacyjnych. Jedynym łączącym strumieniem 
wyjściowym podsystemu jest strumień (częściowo) oczyszczonego w podsys- 
temie gazu. Operator gałęzienia typu R jest więc nieaktywny i przestrzeń 
stanów zawiera 35 stanów. Optymalną strukturę instalacji, wraz ze stężeniami 
składnika zanieczyszczającego w poszczególnych strumieniach, pokazuje ryc. 4. 


200 0.0001 
F, 0.08 1,K 1, ADF —— 
0.02 -  1o02 
| 

|BSAESEE! 

F,=="3H] RZE | 

[R .004 ! 0.004 
F, __0004 rank I 0.00 '3,ABI 

REZ icz 


Ryc. 4. Optymalna struktura instalacji oczyszczania gazów 


Podsystemy zaznaczone linią przerywaną oznaczają podsystemy wyeliminowa- 
ne w wyniku optymalizacji parametrycznej stanu przy obliczaniu wartości 
funkcji oceniającej (13) i nie występują fizycznie w zsyntezowanym systemie. 
Dla porównania, koszt instalacji oczyszczającej mieszaninę strumieni jest 
wyższy o blisko 40% , instalacji zaś, w której każdy strumień oczyszczany jest 
w oddzielnym ciągu — o blisko 30%. 


SYNTEZA WĘZŁA KRYSTALIZACJI 
W INSTALACJI REGENERACJI KĄPIELI KOAGULACYJNEJ 
W WYTWÓRNI WŁÓKIEN WISKOZOWYCH [5] 


W wytwórni włókien wiskozowych powstaje F=390000 kg/h wyczerpanej 
kąpieli koagulacyjnej o składzie: 


XNa1s0, 02  Xpso, 01 Xp,o 0,7 


i temperaturze tp=45'C. 
Z kąpieli tej należy usunąć DEL= 370 kg/h Na,SO4 w postaci soli bezwodnej. 
Należy zsyntezować węzeł krystalizacyjny, w którym jednostkowe zużycie 
energii na 1 kg wydzielonej soli będzie minimalne. 

Na ryc. 5. przedstawiono przykładowy wykres równowagowy układu woda 
(w) — siarczan sodu (SN) — kwas siarkowy (KS) w temperaturze 25C 
i w zakresie stężeń kwasu siarkowego poniżej 70% masowych. Jak widać, 
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H,O XKS HSO, 


Ryc. 5. Przykładowa izoterma równowagowa w układzie H,0 — Na,SO, — HSO, 


w układzie tworzy się pięć soli: 1) dziesięciowodny siarczan sodowy 
Na,SO,:10H,0 (SG), 2) bezwodny siarczan sodowy Na,SO, (SN), 3) sól 
mieszana Na„H(SO,), (SM), 4) uwodniona sól kwaśna NaHSO;, : H,O (SKW), 
5) bezwodna sól kwaśna NaHSO, (SK). 

Etykiety tych soli mogą zostać uporządkowane za pomocą relacji 
zdefiniowanej w następujący sposób: 


isRjs, 


xrKs |;s XSks |js  LSks |js LSKs lis 
dla is, jse (SG, SN, SM, SKW, SK). 


W rozpatrywanym układzie ciało stałe — ciecz mamy: 


SG R SN RSM RSKW RSK. 


jeżeli 


788 480 
pary wodnej pary wodnej 


1872 2568 234 
(0.2,041) (0.34,0.15] (002, 080) 
1066 
NazHISO,)> 
370 
Na SO, 


Ryc. 6. Optymalna struktura węzła krystalizacji 
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Przyjęto, iż węzeł krystalizacji będzie zbudowany z podsystemów krystalizacyj- 
nych, w których będą krystalizowały sole SN, SM i SK. Etykiety podsystemów 
mogą być więc również uporządkowane relacją R. Schemat łączących strumie- 
ni zewnętrznych podsystemu krystalizacyjnego jest identyczny z przedstawio- 
nym na ryc. 2, przy czym indeks r odnosi się do strumienia roztworu, zaś 
indeks s do strumienia oddzielonej fazy stałej. Ze względu na pojedynczą 
surówkę nieaktywny jest operator gałęzienia S i przestrzeń stanów zawiera 36 
stanów. Optymalną strukturę węzła, wraz z niektórymi parametrami, pokaza- 
no na ryc. 6. 


WNIOSKI 


Możliwość uporządkowania etykiet typów podsystemów, które mogą wejść 
w skład syntezowanego systemu, pozwoliła na sformułowanie efektywnego 
algorytmu syntezy opartego na strategii poszukiwawczej w dyskretnej prze- 
strzeni stanów — struktur zupełnych systemu. Działanie algorytmu zostało 
wielokrotnie sprawdzone na przykładach syntezy systemów procesowych 
znacznie różniących się stroną fizyczną, przy różnych liczbach strumieni 
surówkowych i ich parametrach, rozmaitych postaciach funkcji celu, ograni- 
czeniach i uwarunkowaniach problemu. Dodatkową zaletą algorytmu jest 
łatwość wbudowania weń heurystyk odsiewających wynikłych z doświadczenia 
projektanta, wskazówek literaturowych lub analizy wyników syntez przep- 
rowadzonych na układach testowych. Można przypuszczać, że obszar za- 
stosowań algorytmu będzie się nadal rozszerzał. 


"WYKAZ OZNACZEŃ 


C zbiór cech podsystemów 

F _ natężenie przepływu surówki 

G _ zbiór cech strumieni kluczowych 

H _ zbiór binarnych zmiennych strukturalnych kodujących strukturę systemu 
(np. macierz) 

H* zbiór etykiet wyrazów niezerowych zbioru H 

HS zbiór kodujący stan 

IS zbiór etykiet podsystemów wchodzących w skład systemu 

IT zbiór etykiet typów podsystemów 

Nf liczba strumieni surówkowych 

R relacja porządkująca etykiety typów podsystemów; mRn oznacza, że 
podsystem typu m poprzedza bezpośrednio podsystem typu n w uporząd- 
kowanej sekwencji etykiet; mRRn oznacza, że podsystem typu m poprze- 
dza (niekoniecznie bezpośrednio) podsystem typu n 

S$ zbiór cech strumieniowych 

SZ zbiór etykiet struktur zupełnych 

V zbiór zmiennych procesowych 
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x _ ułamek masowy 

xr* skład roztworu nasyconego 

xs skład soli 

y ułamek molowy składnika zanieczyszczającego 


ETYKIETY I INDEKSY 


f sekwencyjny 

F w surówce 

ij etykiety podsystemów; i, je IS 

is,js etykiety soli 

k,l etykiety linii procesowych; k,le(1, 2,..., Nf) 
m,n etykiety typów podsystemów; m,n EIT 

0 etykieta otoczenia systemu 

P,q etykiety strumieni 
8 


roztwór 
re recyrkulacyjny 
s faza stała 
st etykieta typu strumienia; st=r lub st=s 
we wejściowy 
wy wyjściowy 


JAN DUDCZAK 


SYNTHESIS OF THE CHEMICAL ENGINEERING SYSTEMS 
WITH KNOWN RELATION ORDERING THE SUBSYSTEMS 


Summary 


An algorithm of synthesis of the chemical engineering systems whose subsystems can be 
precedence ordered with the help of some ordering relation R has been developed. The algorithm is 
based on the search strategy in a discrete space of states. 

It can be applied to the synthesis of the systems fed with a single feed stream as well as with 
multiple feeds. Two examples of the application of the algorithm to the synthesis of a gas cleaning 
system for multisource emission and of a crystallization system in coagulation bath regeneration in 
manufacture of viscose fibers have been presented. 


AH NYJUAK 


CMHTE3 CHCTEM XMMMUECKOM IIPOMBINJIEHHOCTM 
C M3BECTHBIM COOTHOLIEHHEM YIIOPANOUEHHA 
IIOĄCHCTEM 


Pe3rome 


Pa3pa6oTaH aJlropHTM CHHTE3a CHCTEM XHMH1ECKOi MpOMBIIIIJIEHHOCTH, IMOĄCHCTEMBI 
KOTOPAIX MOTYT OBITb yYIIOPANO4EHBI HPA IOMOIIM ONPEĄEJEHHOTO COOTHOLIEHHA R. ANropuTM 
OCHOBAH Ha CTpATETHK IIOHCKA B NACKPETHOM NpOCIPAHCTBE COCTOAHHIi. OH MOXET HpUMEHATECA 


5 Szczecińskie Rocz Nauk. t VI, z I 


66 


JULA CHHTE3A CHCTEM, IHTAEMBIX OMHHM HJIH HECKOJIEKAMH IIOTOKAMH CHIDbA. IIpe4CTABJIEHbl 
IIpAMepBI IIpAMeHeHHA AJITODHTMA JUIA CHHTE3A YCTAHOBKA O4HCTKA TA30B, HPOHCXOJALNAX H3 
MHOTOHCTOTHHKOBOŃ 3MHCCHH, A TAKOKE CHHTE3A yZJIA KDACTAJLIAZALHH B YCTAHOBKE pereHepau HH 
KOATYJLANHOHHOH BAHHBI JUIA IIDOH3BOĄCTBA BACKOZHBIX BOJIOKOH. 
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JOANNA KARCZ 


DOŚWIADCZALNE BADANIA 
BFEKTYWNOŚCI MIESZANIA CIECZY 
DLA UKŁADU DWU MIESZADEŁ TURBINOWYCH 


Instytut Inżynierii Chemicznej i Chemii Fizycznej Politechniki Szczecińskiej, 
al. Piastów 42, 71-065 Szczecin 


Streszczenie: W artykule przedstawiono wyniki doświadczalnych badań mocy 
mieszania i wnikania ciepła dla ogrzewanego płaszczowo mieszalnika, zaopat- 
rzonego w dwa mieszadła turbinowe na jednym wale. Pomiary przeprowadzono dla 
różnych wartości parametru H/D. Na podstawie kryterium (4) oszacowano efektyw- 
ność mieszania badanych zestawów mieszadeł. Stwierdzono, że z punktu widzenia 
zastosowanego kryterium, najkorzystniejszy dla zbiornika smukłego (H/D > 1) jest 
zestaw dwu mieszadeł turbinowych tarczowych. 


Słowa kluczowe: mieszanie mechaniczne cieczy — moc mieszania — wnikanie 
ciepła — efektywność mieszania — układ dwu mieszadeł. 


WSTĘP 


Mieszanie mechaniczne jest od dawna szeroko stosowane w praktyce 
przemysłowej, jednak w wielu przypadkach brak jest do dziś pełnego opisu 
ilościowego wpływu intensywności mieszania na procesy przebiegające w mie- 
szalniku. Wynika to między innymi z dużej różnorodności stosowanych 
w praktyce konstrukcji tych aparatów. W wielu zastosowaniach technicznych 
wyłania się w związku z tym problem doboru najkorzystniejszego typu 
mieszalnika dla konkretnego procesu mieszania. Kryterium porównawczym 
może być w takim przypadku efektywność mieszania. 

Pojęcie to można formułować bardzo ogólnie jako wielkość świadczącą 
o sprawności wykorzystania energii mieszania dla określonego celu techno- 
logicznego [1]. W najogólniejszym przypadku celem tym może być prze- 
prowadzenie reakcji chemicznej w reaktorze zbiornikowym z mieszadłem 
mechanicznym. Osiągnięcie tego celu może być jednak uwarunkowane przez 
szybkość procesów cząstkowych, zachodzących równocześnie, np. intensyw- 
ność wymiany ciepła czy szybkość homogenizacji. W literaturze znane są więc 
propozycje badania efektywności mieszania przy wykorzystaniu wyżej wymie- 
nionych procesów jako procesów kontrolnych, na podstawie których ocenia się 
jakość pracy danego mieszadła. O tym, które z kryteriów będzie najbardziej 
odpowiednie, decyduje rodzaj reakcji chemicznej. Dla reakcji o dużym efekcie 
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cieplnym kontrolowanie intensywności wymiany ciepła z otoczeniem może 
mieć zasadnicze znaczenie dla oceny pracy danego typu mieszadła. W prezen- 
towanych badaniach przyjęto takie właśnie kryterium do analizy intensywności 
mieszania przy zastosowaniu różnych geometrii układu mieszadło — zbiornik. 

Aby jednak móc obliczyć efektywność mieszania, konieczna jest znajomość 
równań, opisujących moc mieszania i wnikanie ciepła, dla mieszalników 
o różnych parametrach geometrycznych. Ze względu na duże znaczenie 
praktyczne, między innymi zastosowanie w bioinżynierii, szczególnie potrzebne 
wydaje się wypełnienie istniejącej w literaturze luki, dotyczącej danych do- 
świadczalnych dla mieszalnika, zaopatrzonego w więcej niż jedno mieszadło na 
wspólnym wale. 


KRYTERIUM EFEKTYWNOŚCI 


Wyniki badań mocy mieszania i wnikania ciepła przedstawiane są najczęś- 
ciej za pomocą równań kryterialnych [1]. Równanie mocy mieszania, obowią- 
zujące dla burzliwego zakresu przepływu cieczy w mieszalniku, ma następującą 
postać 


P 
AES nalo E (1) 


Wartość współczynnika wnikania ciepła « między ścianką mieszalnika a mie- 
szaną cieczą można określić na podstawie równania 
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w którym wykładniki A=2/3, B=1/3 oraz E=0,14 [1,2]. Występujące 
w równaniach (1) i (2) współczynniki C, są funkcją parametrów geometrycz- 
nych mieszalnika. 

Założona intensywność wnikania ciepła dla mieszalnika o określonej 
geometrii (w naszym przypadku wartość «) jest osiągana kosztem różnej 
wielkości energii mieszania. W tym sensie kryterium intensywności wnikania 
ciepła nie może być wystarczające przy wyborze najkorzystniejszego typu 
konstrukcji mieszalnika. Przyjęte kryterium efektywności powinno mianowicie 
uwzględniać nakład energii, przy którym osiąga się żądaną intensywność 
wnikania ciepła. Warunek ten spełnia funkcja typu 


bi 
1=bo (7) | 0) 


gdzie P/V, jest energią mieszania na jednostkę objętości mieszanej cieczy. 
Postać funkcji (3) powinna być również przeanalizowana z punktu widzenia 

kryteriów powiększania skali mieszalnika, co ma znaczenie dla aparatów 

przemysłowych. Proces wnikania ciepła można bowiem modelować przy 


=G> RE PORY; (2) 
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zachowaniu w obu porównywanych aparatach stałej wartości jednego z warun- 
ków: liczby Re, prędkości obwodowej mieszadła, częstości obrotów mieszadła, 
jednostkowego momentu obrotowego lub jednostkowej energii mieszania. 
Modelowanie według ostatniego z wymienionych kryteriów zapewnia bardzo 
małe zmiany współczynnika wnikania ciepła ze zmianą średnicy mieszalnika 
i z tego powodu jest zalecane w literaturze [3]. Kierując się powyższymi 
względami, w niniejszej pracy zdecydowano przyjąć do analizy własnego 
materiału doświadczalnego kryterium efektywności (3). 


CZĘŚĆ DOŚWIADCZALNA 


Badania mocy mieszania i wnikania ciepła przeprowadzono dla dwu 
standardowych mieszadeł: turbinowego tarczowego (ryc. 1) i turbinowego 
o łopatkach pochylonych (ryc. 2) zamontowanych w mieszalniku z prze- 
grodami, ogrzewanym za pomocą płaszcza. W wyniku wcześniejszej analizy 
doświadczalnej pewne zestawy tych mieszadeł odrzucono, wybierając do 
szczegółowszych badań następujące kombinacje mieszadeł: 


Ryc. 1. Podstawowe parametry geometryczne Ryc. 2. Podstawowe parametry geometryczne 
mieszadła turbinowego tarczowego mieszadła turbinowego o łopatkach pochy- 
lonych 


1) dwa mieszadła turbinowe tarczowe o geometrii standardowej, tzn. 
d/D=0,33; a/d=0,25; b/d=0,2; Z=6, 

2) mieszadło turbinowe tarczowe standardowe (górne) oraz mieszadło 
turbinowe o łopatkach pochylonych (dolne), o następującej geometrii: 
d/D=0,33; b/d=0,2; Z=6; a=135, 

3) dwa mieszadła turbinowe o łopatkach pochylonych (d/D=0,33; b/d=0,2; 
Z=6, «=135 (górne) oraz «=45* (dolne)), 
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4) standardowe mieszadło turbinowe tarczowe (badania przeprowadzono 
w celach porównawczych). 

Pomiary przeprowadzono w wyposażonym w standardowe przegrody mie- 
szalniku o średnicy wewnętrznej D=0,45 m, dla czterech różnych wysokości 
H napełnienia cieczą zbiornika: H=D; H=1,25D; H=1,5$D; oraz H=1,9D. 
Odległość dolnego mieszadła od dna zbiornika wynosiła h, =H/3, a górnego 
h,=2H/3. Mieszano olej maszynowy M6 w zakresie zmian liczby 
Re<10*%; 3:10*>. 

Pomiary użytecznego momentu obrotowego wykonano metodą tensomet- 
ryczną. Badania wnikania ciepła przeprowadzono metodą pomiaru średnich 


13 
12 


Ryc. 3. Schemat aparatury doświadczalnej 


1 — chłodnica oleju, 2 — pompa obiegowa, 3 — zbiornik pośredni, 4 — zbiornik mieszalnika, 5 — przegroda, 6 — mieszadło, 7 — płaszcz, 

8 — manometr, 9 — tuleja skrętna, 20 — silnik elektryczny, /1 — przekładnia pasowa, 12 — magnes, 13 — kontaktrony, 14 — komutator, 

15 — wzmacniacz tensometryczny, 16 — rcjestrator MK6, 17 — zasilacz, 18 — licznik elektroniczny, 19 — termostat, 20 — termoelement, 
2] — urządzenie sterownicze, 22 — garnek kondensacyjny, 23 — rotametr 
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wartości współczynnika wnikania ciepła. Schemat aparatury doświadczalnej 
przedstawiono na ryc. 3. Zbiornik mieszalnika (4) otoczony był płaszczem 
grzejnym (7), który wykonano z blachy miedzianej o grubości 5 mm. 
W mieszalniku istniała możliwość przesuwu dna wzdłuż osi pionowej, co 
pozwalało regulować wysokość H napełnienia cieczą zbiornika. Odległość mie- 
szadeł od dna zbiornika można było zmieniać przez odpowiednie ustawienie 
płyty, na której zamocowano napęd wału głównego. W zbiorniku wyposażo- 
nym w standardowe przegrody (5) montowano odpowiedni zestaw dwu 
mieszadeł (6). Wał główny napędzano za pomocą bocznikowego silnika prądu 
stałego (10). Napęd z silnika na wał przenoszono za pośrednictwem przekładni 
z paskiem klinowym (11). Do regulacji i kontroli clektrycznego układu 
napędowego służyło urządzenie sterownicze (21), wyposażone w regulator 
prądu wzbudzenia, miernik prądu twornika i regulator liczby obrotów silnika. 

Użyteczny moment obrotowy mieszadła mierzono za pomocą torsjometru 
tensometrycznego, zbudowanego we własnym zakresie według zasady zmo- 
dyfikowanej w stosunku do rozwiązania przedstawionego w pracy [4]. 
Głównym elementem torsjometru była wymienna stalowa tuleja skrętna (9) 
o średnicy 30,8 mm oraz o grubości ścianki 0,4 mm. Na powierzchni tej tulei 
naklejano cztery tensometry, połączone w układ mostka. Sygnał z tensometrów 
odprowadzono do wirujących części komutatora (14) (zbieracz prądu typu 
Z12U) przewodami umieszczonymi w wydrążeniu w wale głównym. Sygnał 
z nieruchomych końcówek komutatora przesyłano do wzmacniacza (15) 
(wzmacniacz cyfrowy KWS 3082 A firmy Hottinger). 

Częstość obrotów mieszadła określano za pomocą elektronicznego licznika 
impulsów. Czujnikiem impulsów były kontaktrony (13), zamocowane na 
tekstolitowym pierścieniu, umieszczonym nad wałem przekładni, na którym 
przyklejono magnes (12). Sygnały elektryczne przesyłane były do układu 
współpracującego z elementami stykowymi (17) i wyświetlane w okienku 
licznika elektronicznego (156). 

Krążenie oleju wymuszała pompa obiegowa (2). Do mieszalnika olej 
podawano poprzez zawory regulacyjne, rotametr (23) oraz sitko umieszczone 
w dnie zbiornika. Z przelewów mieszalnika kierowano go do zbiornika 
pośredniego (3), a następnie do przeponowej chłodnicy oleju (1). Olej znaj- 
dujący się w mieszalniku ogrzewano płaszczowo parą wodną. Kondensat 
odprowadzano do garnków kondensacyjnych (22). Pomiaru temperatur doko- 
nywano za pomocą termometrów rtęciowych oraz termoelementów Fe — cons- 
tantan (20) podłączonych do rejestratorów (16). 


WYNIKI BADAŃ I DYSKUSJA 


W wyniku przeprowadzonych badań otrzymano dwie serie pomiarów 
mocy mieszania i współczynnika wnikania ciepła, po 160 punktów doświad- 
czalnych każda. 
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Charakterystyki mocy dla wszystkich badanych układów geometrycznych 
przedstawiono graficznie w postaci funkcji Ne=/f(Re). Na ryc. 4 i 5 pokazano 


| aaa [a + 
[4/0] 1 [12] 15 | 18 | 


250 kożox torożnt xO= 


105 15-104 20-10* 250% Re 


Ryc. 4. Charakterystyka mocy Ne=/ (Re) dla standardowego mieszadła turbinowego tarczowego 
i różnych wartości parametru H/D 


Ne 
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4 404 104 40% 
10 15:10 20:10 25:10 Re 


Ryc. 5. Charakterystyka mocy Ne=/ (Re) dla układu dwu mieszadeł turbinowych o łopatkach 
pochylonych pod kątem «=135? (górne) i «=45” (dolne) i różnych wartości parametru H/D 


wykresy dla mieszalnika zaopatrzonego w jedno standardowe mieszadło 
turbinowe tarczowe oraz dwa mieszadła turbinowe o łopatkach pochylonych 
pod kątem «=45" (dolne) i a =135” (górne). Cechą wspólną tych wykresów jest 
to, że liczba mocy w obu przypadkach nie zależy od parametru H/D, wobec 
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czego w zakresie H/D <1; 1, 9> można przyjąć Ne=const. Na dwu kolejnych 
rysunkach przedstawiono charakterystyki mocy dla następujących zestawów 
mieszadeł: turbinowego tarczowego (górne) i turbinowego o łopatkach po- 
chylonych pod kątem «=135” (dolne) (ryc. 6) oraz dwu mieszadeł turbinowych 
tarczowych (ryc. 7). 


4 4104 404 , 
10 1.510 2010 2540% Re 


Ryc. 6. Charakterystyka mocy Ne=/(Re) dla układu dwu mieszadeł: turbinowego tarczowego 
(górne) i turbinowego o łopatkach pochylonych pod kątem «=135? i różnych wartości parametru 
H/D 


10% 1.510% 20-10% 250% BA 


Ryc. 7. Charakterystyka mocy Ne=/(Re) dla układu dwu standardowych mieszadeł turbinowych 
tarczowych i różnych wartości parametru H/D 
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W obu przypadkach ze wzrostem wartości parametru H/D zwiększają się 
wartości liczby Ne. Wartości liczby mocy dla wszystkich badanych układów 


geometrycznych zestawiono w tabeli 1. 


Tabela 1. Wartości liczby Ne dla badanych zestawów mieszadeł 


H/D 
L Mieszadło Mieszadło 

F: górne dolne 1 1,25 1,5 19 

Ne 
1 — TT 5,81 5,80 5,85 5,96 
2 iE2135> TŁ 45 3,80 3,59 3,45 3,80 
3 TT TŁ 135 6,50 1,09 1,49 17,60 
4 TE TT 8,46 10,61 11,10 11,41 


TT — mieszadło turbinowe tarczowe; TŁ — mieszadło turbinowe o łopatkach pochylonych. 


W przypadku zestawów dwu różnych mieszadeł (poz. 2 i 3 w tab. 1) liczba 
mocy jest w przybliżeniu sumą liczby Ne dla pojedynczych mieszadeł. 
Z porównania liczb mocy dla mieszadła turbinowego tarczowego (poz. 1 w tab. 
1) i układu dwu takich mieszadeł (poz. 4 w tab. 1) wynika, że dla wartości 
parametru H/D=l liczba Ne dla dwu mieszadeł jest mniejsza niż podwojona 
liczba mocy dla pojedynczego mieszadła. W przybliżeniu tę dwukrotną wartość 
liczby Ne dla układu dwu mieszadeł turbinowych tarczowych otrzymuje się 
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2.2 


4.0 41 4,2 43 


lg Re 
Ryc. 8. Zależność Nu/Pr***V"''* =f(Re) dla standardowego mieszadła turbinowego tarczowego 
i różnych wartości parametru H/D 
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4.0 4,1 4,2 4,3 Lg Re 


Ryc. 9. Zależność Nu/Pr*3 V9'1*=f(Re) dla układu dwu mieszadeł turbinowych o łopatkach 
pochylonych pod kątem «=135? (górne) i u=45” i różnych wartości parametru H/D 


przy wartości parametru H/D=1,9. W tych warunkach prawdopodobnie 
mieszadła przestają oddziaływać na siebie hydrodynamicznie. 

Rezultaty badań współczynnika wnikania ciepła przedstawiono w postaci 
funkcji Nu/Pr033 Vo'1*=f(Re) na ryc. 8-11. Analiza wykresów na ryc. 
8 i 9 wskazuje, że w przypadku mieszalnika zaopatrzonego w jedno mieszadło 
turbinowe tarczowe lub dwa mieszadła turbinowe o łopatkach pochylonych 


2,4 


2.2 


1.0 41 4.2 4.3 lg Re 


Ryc. 10. Zależność Nu/Pr33 V":!*=/(Re) dla układu dwu mieszadeł: turbinowego tarczowego 
(górne) i turbinowego o łopatkach pochylonych pod kątem «= 135? i różnych wartości parametru 
H/D 
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największa intensywność wnikania ciepła występuje dla wartości parametru 
H/D=1 i bardzo wyraźnie maleje ze wzrostem parametru H/D. Jak wynika 
natomiast z ryc. 10, wzrost wartości parametru H/D powoduje niewielkie 
zmniejszenie intensywności wnikania ciepła dla mieszalnika zaopatrzonego 
w dwa mieszadła: turbinowe tarczowe (górne) oraz turbinowe o łopatkach 
pochylonych pod kątem «=135” (dolne). Zamontowanie dwu mieszadeł 
turbinowych tarczowych na jednym wale powoduje wzrost intensywności 
wnikania ciepła w miarę wzrostu wartości parametru H/D, co potwierdzają 
wykresy, pokazane na ryc. 11. 
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40 4 4.2 4,3 tg Re 


Ryc. 11. Zależność Nu/Pr'*33 V"'!* =f/(Re) dla układu dwu standardowych mieszadeł turbinowych 
tarczowych i różnych wartości parametru H/D 


Tabela 2. Wartości stałych C, w równaniu (2) dla badanych zestawów mieszadeł 


H/D 
Mieszadło Mieszadło 
Lp. górne dolne 1 1,25 1,5 1,9 
C, 
1 — FT 0,732 0,678 0,601 0,547 
2 TŁ/1357 TŁ 45 0,731 0,657 0,562 0516 
3 TT TŁ 1359 0,738 0,729 0,690 0,661 
4 m TE 0,702 0,842 0,784 0,753 


Wyniki pomiarów współczynnika wnikania ciepła opracowano w posta- 
ci równania kryterialnego (2), a wartości współczynnika C, zostały zebrane 
w tab. 2. 
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Na podstawie danych doświadczalnych obliczono również dla każdej 
badanej geometrii mieszalnika wartości stałej by oraz wykładnika b,w rów- 
naniu (3), przyjętym jako kryterium efektywności. Stwierdzono, że w przedziale 
ufności obliczonych wartości wykładników b, mieści się wartość b, =0,22 
i w związku z tym tę właśnie wartość wspólną dla wszystkich badanych 
geometrii mieszalników przyjęto do dalszych obliczeń. Wyniki badań aprok- 
symowano wobec tego następującą funkcją 


0,22 
EX Q | (4) 


Wartości stałej by dla badanych geometrii mieszalników zebrano w tab. 3. 
Porównanie wartości stałych b, prowadzi do interesujących wniosków. Wśród 
porównywanych geometrii mieszalników, charakteryzujących się standardową 
wartością parametru H/D (H/D=1) najkorzystniejszy w aspekcie przyjętego 
kryterium okazuje się zestaw dwu mieszadeł turbinowych o łopatkach po- 
chylonych pod kątem «a=135” (górne) i u =45" (dolne). 


Tabela 3. Wartości stałej by w równaniu (4) dla badanych zestawów mieszadeł 


H/D 
Mieszadło Mieszadło 
Lp. górne dolne 1 1,25 1,5 19 
bą 
1 — TT 168,90 164,66 150,63 146,82 
2 TŁ 1353 TŁ 45 185,77 176,87 158,82 152,24 
3 TT TŁ "1355 166,88 169,68 166,49 168,69* 
4 TE TT 147,28 179,52 172,77 175,64 


* Wartość średnia by= 167,94 dla H/D<1; 19>. 


Dla zestawu mieszadeł: turbinowe tarczowe (górne) i turbinowe o łopatkach 
pochylonych pod kątem «a =135” (dolne) wartość stałej b, nie zmienia się 
w całym zakresie zmian parametru H/D. 

W przypadku zbiornika o dużej smukłości najlepsze rezultaty osiąga się 
stosując zestaw dwu mieszadeł turbinowych tarczowych. 

Otrzymane na podstawie danych doświadczalnych wykresy funkcji 
«=f(P/V,) przedstawiono również na ryc. 12-15. Linie ciągłe na tych rysun- 
kach odpowiadają równaniu (4). Przebieg funkcji na ryc. 12-15 jest graficznym 
potwierdzeniem wniosków, wynikających z tabeli 3. 

Reasumując, przeprowadzone obliczenia efektywności mieszania wskazują 
na to, że różne układy geometryczne mieszadeł wykazują znacznie różną 
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e. 50458. 7-884 2 3 4 5 6407 Pfw 
w Lr] 
Ryc. 12. Zależność u =f/(P/V,) dla mieszalnika zaopatrzonego w standardowe mieszadło turbinowe 
tarczowe 


Ryc. 13. Zależność a=/(P;V,) dla mieszalnika zaopatrzonego w dwa mieszadła turbiiowe 
o łopatkych pochylonych pod kątem «=135* (górne) i u=45" (dolne) 


efektywność. Jest to związane z charakterem cyrkulacji cieczy, która jest 
wytwarzana przez te zestawy mieszadeł. Przy założonej energii mieszania 
najkorzystniejsze warunki wymiany ciepła zapewnia taki zestaw mieszadeł, za 
pomocą którego osiągane są największe prędkości w pobliżu powierzchni 
wymiany, jak również najwyższy poziom burzliwośći w obszarze przyściennym 
mieszalnika. 


-—-.1—| 


3. 4 5 6 7 8 9 102 2 3 % 5406 ABT 
[ns] 


c 
Ryc. 14. Zależność «= (P/V,) dla mieszalnika zaopatrzonego w dwa mieszadła: turbinowe tarczowe 
(górne) i turbinowe o łopatkach pochylonych pod kątem «= 135" (dolne) 


2 | 


3 4 5 6 7 8 9102 2 4 5407 v [5] 


Ryc. 15. Zależność u=/(P/V.) dla mieszalnika zaopatrzonego w dwa standardowe mieszadła 
turbinowe tarczowe 


PODSUMOWANIE 


1. Na podstawie doświadczalnych badań mocy określono wartości liczby 
Ne dla czterech wybranych zestawów mieszadeł w zakresie burzliwego prze- 
pływu cieczy w mieszalniku. 

2. Wyniki doświadczalnych badań współczynnika wnikania ciepła uogól- 
niono w postaci równania kryterialnego (2). 

3. Wykazano doświadczalnie, że dla standardowej wartości parametru 
H/D=l zestaw dwu mieszadeł turbinowych o łopatkach pochylonych pod 
kątem «=45” (dolne) oraz «=135” (górne) zapewnia lepszą efektywność 
mieszania cieczy niż pojedyncze mieszadło turbinowe tarczowe. 
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4. Stwierdzono, że w aspekcie przyjętego kryterium efektywności najkorzy- 
stniejszy dla smukłych zbiorników (H/D> 1) jest układ dwu mieszadeł tur- 
binowych tarczowych. 


WYKAZ OZNACZEŃ 


A,B,E [-] wykładniki w równaniu (2) 

a [m] długość łopatki mieszadła 

b [m] szerokość łopatki mieszadła 

bo [-] stała w równaniu (3) 

by [-] wykładnik w równaniu (3) 

Cys 6; [-] współczynniki w równaniach (2) i (3) 

c [Jskg''-K"!] ciepło właściwe mieszanej cieczy 

D [m] wewnętrzna średnica mieszalnika 

d [m] średnica mieszadła 

H [m] wysokość słupa cieczy w mieszalniku 

h, [m] wysokość zawieszenia mieszadła liczona od dna 
mieszalnika do osi symetrii mieszadła 

n [1/s] częstość obrotów mieszadła 

P [W] moc mieszania 

V, [m*] objętość mieszanej cieczy 

Z [-] liczba łopatek mieszadła 

a [] kąt pochylenia łopatki mieszadła 

a [W:'m*:K"!'] współczynnik wnikania ciepła 

q [kg'm'1:s'!] dynamiczny współczynnik lepkości cieczy 

A [W:'m"':K"!] współczynnik przewodzenia ciepła cieczy 

p [kg/m?] gęstość cieczy 


MODUŁY BEZWYMIAROWE 
Nu=aD/a liczba Nusselta 


Pr=cy/a liczba Prandtla 
Re=nd?p/q liczba Reynoldsa 
V=n/q, simpleks lepkości 


Ne=P(n*d”p) liczba mocy 


JOANNA KARCZ 


EXPERIMENTAL STUDIES 
OF LIQUID MIXING EFFECTIVENESS 
FOR DUAL SYSTEMS OF TURBINE IMPELLERS 


Summary 


The results of the experimental investigations of the power consumption and heat transfer for 
jacketed heated stirred tank equipped with a dual system of turbine impellers are presented in the 
paper. The measurements were carried out for various values of the geometrical parameter H/D. 
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On the basis of the criterion (4) the mixing effectiveness was evaluated for tested systems of 
impellers. 

From point of view of the used criterion, it has been proved that the dual system of two disc 
turbines is the best for slender tank (H/D> 1). 


MOAHHA KAPU 


OKCIIEPAMEHTAJIBHBIE KCCJIENOBAHMA 
ODOEKTHBHOCTH IIEPEMENMABAHMA KHJIKOCTA 
NJIA CAHCTEMPBI Ą|BYX TYPEHHHDIX MEIIAJIOK 


Pe3rome 


B crarbe npeqcTaBJIeHkI pe3yJIbTATKI 3KCIEDHMEHTAJILHAIX HCCJIE1OBAHHHi MOLNHOCTH 
repeMelIABAHHA AH TENJIOOTĄAUH MIA CMECHTEJA C pyGalikoń o6orpeBa, cHaOxeHHOro MByYMA 
TYpÓHAHBIMAH MELIAJIKAMA Ha ONHOM BaJlry. HI3MepeHHa BBIIIOJIHCHBI JULA PAZHKIX 3HAdeHH 
napamerpa H/D. Ha ocHoBaanau KkpuTepna (4) nponsBeqeHa OleHka 3QdhekTABHOCTH 
IepeM€LIHBAHHA JUIA ACCJIENYEMBIX KOMILIEKTOB MELIIAJIOK. 

YVCTAaROBJIEHO, € TOYKA 3DEHHA NIDAHATOTO KPATEDHA, 4TO CAMO BEITONHOŃ JUIA CMECHTEJIA 
H/D>l aBnaeTca KOMIUIEKT, COCTOAINAŃ H3 NBYX TYDÓHHHBIX HACKOBKIX MEIIIAJIOK. 
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DYSKRETNA POLIOPTYMALIZACJA TYPOSZEREGU 
PRĘTÓW KRATOWNICY PRZESTRZENNEJ 


Instytut Inżynierii Lądowej Politechniki Szczecińskiej, 
al. Piastów 50, 70-311 Szczecin 


Streszczenie. W pracy sformułowano i rozwiązano dyskretne zadanie poliop- 
tymalizacji ortogonalnej kratownicy przestrzennej podpartej punktowo w narożach. 
Przeanalizowano dwie zmienne decyzyjne: typoszereg przekrojów poprzecznych 
prętów oraz liczbę elementów wykreślanych z typoszeregu podstawowego. Zbudo- 
wano 6 typoszeregów przekrojów wykorzystując katalog walcowanych rur stalo- 
wych. Liczebność typoszeregów zmieniano w granicach od 4 do 25. Jako kryteria 
optymalizacji przyjęto masę prętów kratownicy, przemieszczenie środkowego węzła 
oraz liczebność typoszeregu prętów. Zbiór rozwiązań kompromisowych zadania 
polioptymalizacji wyznaczono metodą systematycznego przeszukiwania dyskretnego 
obszaru dopuszczalnego. Rozwiązanie preferowane określono metodą kryterium 
globalnego. 

W pracy określono zasady budowy typoszeregu przekrojów o ograniczonej 
liczebności przydatnego w optymalnym kształtowaniu kratownic przestrzennych. 


Słowa kluczowe: polioptymalizacja — kompromis — preferowany — kratow- 
nica — typoszereg — pręt — przekrój — dyskretny. 


WPROWADZENIE 


W zadaniach optymalizacji konstrukcji wg kryterium minimum objętości 
materiału jako zmienne decyzyjne przyjmowane są między innymi wielkości 
opisujące geometrię przekroju poprzecznego prętów. Są to wartości pól 
przekrojów prętów [1—4] lub zmienne opisujące wymiary przekroju [5— 7]. 
Niekiedy zmienną decyzyjną jest liczebność typoszeregu przekrojów [8, 9] lub 
typ przekroju poprzecznego prętów [8— 10]. Ze względu na charakter zmien- 
nych decyzyjnych wyróżnia się optymalizację ciągłą i dyskretną [11]. Op- 
tymalizacja ciągła, najchętniej stosowana w pracach teoretycznych [3, 12], 
umożliwia określenie wymiarów przekrojów ze zbioru liczb rzeczywistych. 
W przypadku konieczności wskazania propozycji projektu technicznego kon- 
strukcji,i w końcowej fazie projektowania przyjmuje się jako optymalne 
wartości dyskretnych zmiennych decyzyjnych najbliższe rozwiązaniom ciąg- 
łym. W wielu rozwiązaniach jest to metoda dająca poprawne rezultaty [4, 13]. 
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Ogólnie jednak wynik optymalizacji ciągłej nie musi być zgodny z wynikiem 
optymalizacji prowadzonej w zbiorach skończonych [1, 14]. 

Optymalizacja dyskretna wymaga określenia explicite wartości zmiennych 
decyzyjnych. Od sposobu ich przyjęcia zależy wynik zadania [8, 15, 16]. 
Wyróżnić można, przy tym dwa rodzaje zmiennych decyzyjnych — związane 
i swobodne. Pierwsze mogą przyjmować tylko z góry określone wartości. Do 
zmiennych związanych należą właściwości mechaniczne materiału, sposób 
podparcia konstrukcji, liczba warstw kratownicy, gęstość siatki prętów, 
konfiguracja prętów, system konstrukcyjno-montażowy, typ przekroju poprze- 
cznego prętów itp. Dyskretne zmienne swobodne mogą być przyjęte na wiele 
różnych sposobów, w zależności od aktualnych potrzeb i możliwości projek- 
towych i wykonawczych. Należą tu między innymi takie zmienne jak wysokość 
konstrukcyjna kratownicy, wielkość pola przekroju poprzecznego prętów, 
zakres jednakowych stref sztywności prętów, wielkości opisujące przekrój 
prętów i węzłów, typoszereg węzłów i jego liczebność. W niniejszej pracy 
przeanalizowano dwie ostatnie z wymienionych zmiennych swobodnych w za- 
stosowaniu do ortogonalnej kratownicy przestrzennej. 

Optymalne kształtowanie kratownicy przy określonym typie przekroju 
poprzecznego prętów może być dokonane trzema sposobami. Pierwszy polega 
na operowaniu pełnym asortymentem profili, co np. dla rur walcowanych daje 
398 elementów, dla kątowników 60 elementów, a dla kwadratowych rur 
ciągnionych na zimno 49 elementów, [17]. Liczba elementów gwałtownie 
wzrasta dla przekrojów złożonych. W przypadku dwóch kątowników równo- 
ramiennych wynosi 1770. Tworzenie takich licznych typoszeregów prętów, 
mimo możliwości zautomatyzowania niektórych czynności (np. progra- 
my JC 12 i JC 12R [18]) jest pracochłonne i w efekcie powoduje nadmierne 
wydłużenie czasu projektowania konstrukcji. 

Drugi sposób polega na wyborze typoszeregu prętów o ograniczonej 
liczebności. Należy przy tym wziąć pod uwagę dwa przeciwstawne czynniki. 
Z jednej strony asortyment powinien być możliwie bogaty, aby wykorzystanie 
poszczególnych przekrojów było maksymalne. Z drugiej jednak strony liczba 
profili powinna być niewielka, aby koszt wykonania konstrukcji był jak 
najmniejszy. W różnych pracach można spotkać liczebność typoszeregu prętów 
od 3 do 33 [2,8, 16, 18]. 

Trzeci sposób polega na projektowaniu nowych wymiarów przekrojów 
dających minimum zużycia materiału przy zachowaniu nośności i warunków 
technologii wytwarzania jak dla profili aktualnie produkowanych [11]. Tego 
typu podejście stosowane jest przy tworzeniu programu walcowania nowych 
profili lub ulepszaniu istniejących. 

W dalszym ciągu rozpatrywane będzie drugie stanowisko, tj. optymalizacja 
kratownicy przy zmieniającym się typoszeregu prętów o zmiennej liczebności. 
Maksymalna liczba elementów w typoszeregu wynosi 25, minimalna 4. Celem 
pracy jest poszukiwanie zasad budowy typoszeregu przydatnego w optymaliza- 
cji kratownic przestrzennych. 
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PRZEDMIOT OPTYMALIZACJI 


W pracy przeprowadzono analizę optymalizacyjną i polioptymalizacyjną 
ortogonalnej kratownicy przestrzennej pokazanej na ryc. 1. Konstrukcja 
oparta jest w narożach przegubowo na słupach o rozstawie 24 x 24 m. Pręty 
kratownicy projektowane są z rur walcowych na gorąco [17] ze stali R35. 
Odległości międzywęzłową a=3,0 m i międzywarstwową h=2,1 m przyjęto 
w wyniku rozwiązania zadania polioptymalizacji kratownicy [19]. Konstruk- 
cja składa się z 512 prętów połączonych przegubowo w 145 węzłach. 


UKŁAD PRĘTÓW | WARSTWA GÓRNA KZĄ 
: LI 1% 18. 13. HH 
WZOWZZCYA 
2SPSPKPKPa, Ba Ba Ba P 
/ 
me gk 
zł 
NOTE DRA 5 EL ""rA 
ż L/2 « 42000 i L/2 = 42000 Ł h= 2400, H 


Ryc. 1. Układ prętów kratownicy przestrzennej z rozkładem sztywności dla rozwiązania preferowa- 
nego 


Na kratownicę działa łączne obciążenie pionowe 1926 N/m?. Składowymi 
obciążenia są: masa kratownicy 226 N/m?, masa pokrycia dachowego 500 
N/m?, obciążenie śniegiem 720 N/m? i obciążenie technologiczne 480 N/m?. 
Założono, że obciążenia od pokrycia dachowego i od śniegu przekazywane są 
w postaci sił skupionych na węzły warstwy górnej. Masę kratownicy oraz 
obciążenie technologiczne rozdzielono na węzły warstwy górnej i dolnej. 
Charakterystyczne wielkości obciążeń przypadające na węzły środkowe oraz 
przyjęte wartości współczynników obciążeń zestawiono na ryc. 2. 


K P,= 1.4 x 6480 N 
R P= 4.2 * 4500 N 
s. 


P, = 44* 2460 N 


R 54.4 * 4047 N 


L/2 = 42000 


j0—0 


pozę 
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Ryc. 3. Analizowana część układu z warunkami brzegowymi 


Ze względu na symetrię konstrukcji i obciążeń analizowano 1/4 część kratow- 
nicy. Przyjęte warunki brzegowe na osiach symetrii pokazano na ryc. 3. 
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SFORMUŁOWANIE PROBLEMU 


Analizę zagadnienia optymalizacji przeprowadzono w dwóch etapach. 
W pierwszym poszukiwano typoszeregu przekrojów pierścieniowych o ograni- 
czonej liczebności (J1 =25), dla którego masa prętów kratownicy przypadająca 
na 1 m? rzutu poziomego, zwana dalej wskaźnikiem zużycia materiału, osiąga 
minimum. W drugim etapie rozwiązano dwa zadania polioptymalizacji, 
przyjmując jako składowe wektorowej funkcji celu wskaźnik zużycia materiału, 
liczebność typoszeregu prętów oraz maksymalne przemieszczenie węzła krato- 
wnicy. Poszukiwano zbioru rozwiązań kompromisowych minimalizującego 
analizowane funkcje. 

Zadanie optymalizacji można sformułować następująco: 


510 
y (f)= min y(T)= min (z », 141) (1) 
TeXo TeXo i=1 
gdzie 
y(T) — wskaźnik zużycia materiału w kg/m”, 
o  — masa właściwa stali, yy =7850 kg/m”, 
A; — pole przekroju poprzecznego i-tego pręta kratownicy, 4, € T, 
l,  — długość i-tego pręta kratownicy, 
L  — rozpiętość kratownicy, I=24 m, 
T  - typoszereg przekrojów minimalizujący funkcję celu, 
T  - wektor typoszeregów przekrojów, TES, 
S  —- katalog rur stalowych bez szwu walcowanych na gorąco [17], 
X, — obszar dopuszczalny zadania wyznaczony przez ograniczenia zestawio- 
ne w tab. 1. 
W tab. 1 przyjęto następujące oznaczenia: 
o,, — naprężenia w i-tym pręcie ściskanym, 
o,  — naprężenia w i-tym pręcie rozciąganym, 
N, — siła w i-tym pręcie kratownicy, 
R  - wytrzymałość obliczeniowa stali R35, 
A,  — smukłość i-tego pręta, 


i,  — promień bezwładności i-tego pręta, 
m,„; — współczynnik wyboczeniowy i-tego pręta wg PN-80/B-03200, 
[K], — macierz sztywności kratownicy w n-tej interacji, 


[0] — wektor przemieszczeń węzłów, 

[P] — wektor obciążeń węzłów, 

Owa, — Maksymalne przemieszczenie pionowe węzłów kratownicy, 
— liczba iteracji projektowania, 

M  - liczba analizowanych typoszeregów prętów, 

g — grubość ścianki pręta, 

D  — średnia zewnętrzna przekroju pręta, 

J . — liczebność typoszeregu przekrojów, 


g,  — katalog grubości [17], 
D, — katalog średnic zewnętrznych [17]. 
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Tabela 1. Ograniczenia obszaru dopuszczalnego X, 


Ograniczenia projektowe Ograniczenia technologiczne 
Nim, 2,9XKq<8,8 mm (9) 
0, = <R=210 MPa (2) 
i 38,0 < D<193,7 mm (10) 
N, 3,24 A, 440,3 cm? (11) 
0,, =— <R=210 MPa (3) 
A, DZ ET Me: TM. TEM] (12) 
l Tm=fAy, ..., A,,..., A 13 
1,=—<250 (4) (A, j „) (13) 
, jef=f1,2,..,J) (14) 
[K]„[6]=[P] 6) | J=(25) (15) 
(BŁ g € g.=(2,9; 3,2; 3,6; 4,0; 4,5; 5,0; 
Baz 5 PLM (6) 5,6; 6,3; 7,1; 8,0; 8,8; ...; 
30,0) (16) 
Ograniczenia obliczeniowe D e D,= 38,1; 33,7; 38,0; 42,4; 44,5; 
48,3; 51,0; 54,0; 57,0; 60,3; 
N<10,neW=f1,2,..,N) (7) 63,5; 70,0; 76,1; 88,9; 101,6; 
108,0; 114,3; 133,0; 139,7; 
u A 159,0; 168,3; 193,7; ...; 508,0) (17) 
TS. (18) 


Ograniczenia (12) -(18) powodują, że obszar dopuszczalny zadania X, jest 
zbiorem dyskretnym. Dopuszczalne wartości g i D pokazano na ryc. 4. 
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Ryc. 4. Dopuszczalne wartości grubości ścianek g i średnic zewnętrznych 
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W zadaniach polioptymalizacji analizuje się dwie wektorowe funkcje celu 
JL, (x)=[v(2,), J(x,)] (19) 
f2, (x,)=[y (7), ómsz(€7)]- (20) 


Składowymi wektoróww funkcji celu są: wskaźnik zużycia materiału y(x,), 
liczebność typoszeregu prętów J (x„) oraz przemieszczenie środkowego węzła 
warstwy górnej kratownicy 6,,,(%,). 

Analizuje się dwie zmienne decyzyjne: typoszereg prętów T oraz liczbę 
elementów k wykreślanych z typoszeregu podstawowego. Sposób budowy 
typoszeregów podstawowych i pochodnych omówiono w dalszej części pracy. 
Wektor zmiennych decyzyjnych ma postać 


x,=[T,k]. (21) 


Obszar dopuszczalny w zadaniach polioptymalizacji jest rozszerzony w sto- 
sunku do obszaru X,(X9c X) przez zmianę ograniczenia (15), które przyjmuje 
postać 

J=(25,1357,4b (22) 


Zmianie ulega również liczba analizowanych typoszeregów z M=5 na M=6 
(ograniczenie 8). 


ZBIÓR TYPOSZEREGÓW PRĘTÓW 


Zbiór typoszeregów prętów ustalono po przeprowadzeniu wstępnego 
projektowania konstrukcji. W analizie zastosowano programy DAGI1 i JC12 
[18] wykorzystujące metodę przemieszczeń. Przekroje prętów dobierane są 
automatycznie przez program JC12 z zadeklarowanego typoszeregu. Proces 
projektowania przebiega iteracyjnie. W każdej iteracji dla określonej macierzy 
sztywności [K], wyznacza się siły w prętach kratownicy i wybiera minimalne 
przekroje spełniające warunki (2)—(5) i (13). Zgodność przekrojów w dwóch 
kolejnych iteracjach lub przekroczenie zadeklarowanej liczby iteracji N kończy 
proces projektowania prętów i wyznaczane są przemieszczenia węzłów dla 
obciążeń obliczeniowych. 

W projektowaniu wstępnym przyjęto typoszereg TS składający się z 39 
przekrojów od g,/D, =2,9/31,8 mm, A, =2,63 cm? do g39/Dz,=8,0/219,1 mm, 
Ajzy=53,1 cm?. Obliczenia przerwano po N=$5 iteracjach. Ekstremalne siły 
w prętach wynosiły + 661,8 kN i —477,3 kN. Wymagało to zastosowania 
przekrojów dla rozciągania 6,3/168,3 mm, 32,1 cm? oraz dla ściskania 5,6/159,0 
mm, 27 cm?. Warunek (4) spowodował, że jako minimalny wybrany został 
przekrój 3,2/38,0 mm, 3,50 cm?, dla którego smukłość pręta o długości 3,0 m 
wynosi 1=242,8. Dla projektu wstępnego osiągnięto wartości funkcji celu 
y(TS)=17,319 kg/m?. 

W wyniku przeprowadzonej analizy wstępnej określono minimalny i mak- 
symalny przekrój typoszeregu oraz skorygowano obciążenie masą własną 
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Zmienność pola powierzchni analizowanych typoszeregów 


Ryc. 5. 
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Ryc. 6. Wartości g i D dla typoszere 
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kratownicy. Uwzględniając warunek (4) przyjęto g,/D, =2,9/38,0 mm, A, = 3,20 
cm?, A= 240,9. Dla zabezpieczenia prawidłowego działania programu JC12 
przyjęto gąs/D,s =8,0/168,3 mm, 4,5=40,3 cm?, A„„= 52,9. Liczebność typo- 
szeregów J=25 dobrano tak, aby łatwo można było tworzyć typoszeregi 
pochodne przez pozostawienie co (k+1)-ego przekroju przy założeniu, że 
pierwszy i ostatni zawsze pozostają. W tab. 2 pokazano jak z typoszeregu 
podstawowego (k=0) tworzy się typoszeregi pochodne o zmniejszonej liczebno- 
ści J. Numery przekrojów j w typoszeregach pochodnych wpisane zostały pod 
numerami typoszeregu podstawowego odpowiadającymi jednakowej wielkości 
przekroju poprzecznego A, W zadaniu analizuje się część typoszeregów 
pochodnych pokazanych w tab. 2, spełniających ograniczenie (22). 


Tabela 2. Tworzenie typoszeregów pochodnych z typoszeregu podstawowego 


k je/=(1,...,J) J 
o |1|2|3|4|5|6|7|8|9 |olit|12|13|14|15|16|17|18|19|20|21|22|23|24| 25 
my i 31 14] 151 LólI7] [BI je] (bi |] 2 13 
eAa ż 3 4 5 6 7 8 9 
3 |1 2 3 4 5 6 7 
SA 3 3 4 5 
17 |1 2 3 4 
Trzj| 2 3 
23 |1 ż 


Przy tworzeniu zbioru typoszeregów prętów przyjęto następujące kryteria: 
— każdy typoszereg składa się z J=25 elementów, 
— minimalny i maksymalny przekrój jest identyczny dla wszystkich typo- 
szeregów, 
— wymiary g i D przekrojów należących do typoszeregu spełniają ograniczenia 
technologiczne (9), (10), (16) i (17), 
— wybierane są tylko rury bez szwu walcowane na gorąco, 
— typoszeregi uporządkowane są wg wzrastającego pola powierzchni, przy 
czym zmienia się ono zgodnie ze z góry określoną krzywą (ryc. 5). 

W pracy zbudowano M=$5 typoszeregów przekrojów pierścieniowych. 
Przekroje typoszeregów T1 i T2 dobrane są zgodnie z formułą 


Aj = 1,1113 Aj. (23) 


W typoszeregu T2 utworzono pięć grup przekrojów o jednakowej średnicy 
zewnętrznej i zmieniającej się grubości ścianki (ryc. 6). Pola przekrojów 
typoszeregu T3 wzrastają proporcjonalnie w całym zakresie zmienności. 


Aj 1 =4;+ 1,5458 cm?. (24) 
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Na zasadzie proporcjonalności wzrostu pola przekroju zbudowane są typosze- 
regi T4i T5. W przypadku T4 określono dwa przedziały proporcjonalności 


Aj41=4;,+0,8833 cm? dla j<20, 
Aj41 4A;+4,1750 cm? dla j>20, 

a w przypadku typoszeregu T$ trzy przedziały proporcjonalności 

- Ają124;+2,6750 cm? dla j<5, 

Aj,1 4,+0,9588 cm? dla 5<j<22, (26) 
Ajj1=4,+3,3667 cm? dla j>22. 

Wartości liczbowe średnicy zewnętrznej D i grubości ścianek g dla przyjętych 

typoszeregów pokazano na ryc. 6. 


(25) 


ANALIZA OPTYMALIZACYJNA 


Projektowanie kratownic przeprowadzono programem JC12 na mikrokom- 
puterze typu IBM PC/XT z pamięcią RAM rozszerzoną do 1MB. W danych 
zadeklarowano prowadzenie obliczeń do uzyskania zgodności przekrojów 
w dwóch kolejnych iteracjach lub przerwanie obliczeń po N=10 iteracjach. 
Liczba iteracji uzależniona od trafności wyboru przekrojów startowych wynio- 
sła kolejno dla typoszeregu T1I-N=2, T2—N=10, T3-—N=4, T5S—N=3. 


a) b) 


| Ngs [kN] dla L= 2.25 m a z: Ngs[kN] dla. L= 3.0m ta 
cm 


A 
[re] 


RZA RENE ZEE AREN ME SZ RE 


4 3 5 7 3 4 43 ;, 4 3 5 7 3 44 43 % 


Ryc. 7. Wartości wskaźnika zużycia materiału i przemieszczenia środkowego węzła kratownicy 
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Wyniki analizy przedstawiono na ryc. 7. Na wykresie zaznaczono wartości 
funkcji celu y(T) oraz przemieszczenia środkowego węzła warstwy górnej 
kratownic 0,,,' Minimum masy prętów uzyskano dla typoszeregu T1 

y(£)= min y(T)=y(T,)=16,748 kg/m?. (27) 
p Te Xo 
W typoszeregu T1 każdy następny pręt ma pole przekroju większe od 
poprzedniego blisko 11%. Właściwość ta powoduje dobre wykorzystanie 
naprężeń we wszystkich prętach kratownicy. Musi tu być jednak spełniony 
dodatkowy warunek. Przekroje typoszeregu powinny być tak dobrane, aby 
wzrostowi rozciągającej siły granicznej 


towarzyszył jednocześnie równomierny wzrost ściskanej siły krytycznej 
RA, 


wj 

Ma to decydujące znaczenie przy projektowaniu prętów ściskanych. Procedura 
projektująca zawarta w programie JC12 przeszukuje typoszereg prętów od 
początku, wybierając pierwszy przekrój, dla którego wartość siły krytycznej 
jest większa od obliczonej. Przekroje nie spełniające warunku jednoczesnego 
wzrostu siły granicznej i krytycznej mogą być wybierane dla prętów roz- 
ciąganych, nie są natomiast nigdy wybierane dla prętów ściskanych. 


+ (T) | 6 max 


TS T4 T2 T3 TH T5 M M2 T6 


Ryc. 8. Zależności pola przekroju, siły krytycznej i smukłości prętów 


Przykład tak utworzonego typoszeregu [20] pokazano na ryc. 8a. Prze- 
kroje 4,8,9 są dobrane niewłaściwie dla prętów Ściskanych. Wzrost pola 
przekroju nie jest tu wystarczająco duży w porównaniu ze wzrostem smukłości 
pręta i siła krytyczna maleje zamiast wzrastać. 
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Warunku równomiernego wzrostu wartości siły granicznej i krytycznej nie 
spełnia typoszereg T2 (ryc. 8b). Dla przekrojów od 1 do 5 siła krytyczna 
wzrasta średnio o 9% , a między przekrojami 5 i 6 ponad 5-krotnie. W efekcie 
przekroje 2, 3,4 i 5 ani razu nie zostały wybrane dla prętów ściskanych. 

Najgorszym w sensie minimalizacji wskaźnika zużycia materiału okazał się 
typoszereg T3 o proporcjonalnym wzroście pola powierzchni (ryc. 5). W po- 
równaniu z typoszeregiem T1 wartość funkcji celu zwiększyła się o 10,7%. 

Opracowane typoszeregi prętów porównano z typoszeregiem stosowanym 
przez znaną firmę MERO. Oryginalny typoszereg MERO [21], oznaczony na 
ryc. 5i7 jako M1 daje minimum funkcji celu, ale nie spełnia ograniczeń (9)-(11), 
(15)-(17) i (19) (ryc. 9). Najczęściej wybieranym w projektowaniu przekrojem 


9 $ [mm] 


Ryc. 9. Wartości g i D dla typoszeregów Ml i 16 


jest pierwszy spełniający warunek (4). Z typoszeregu M1 przekrój j=2 nie 
spełniający ograniczeń (9), (16) -(18) jest wybrany dla 188 prętów kratownicy. 
Typoszereg M2 zbudowano z typoszeregu M1 przez wyeliminowanie prze- 
krojów j=l i j=2. W wyniku projektowania kratownicy uzyskano (ryc. 7) 

y(M2)=17,143 kg/m”>y(T1)=16,748 kg/m?. (30) 


Wykorzystując doświadczenia nabyte w trakcie formułowania i analizy 
przedstawionego problemu zbudowano typoszereg T6 spełniający przyjęte 
ograniczenia technologiczne (9)-(18). Wartości D i g dla poszczególnych 
elementów typoszeregu zaznaczono na ryc. 9, a pola przekrojów na ryc. 5. 
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Typoszereg T6 daje najmniejszą, pomijając M1, wartość wskaźnika zużycia 
materiału (ryc. 7). W porównaniu z najlepszym typoszeregiem T1 ze zbioru 
T=(T1,..., T5) wartość funkcji celu (1) zmniejszyła się o 1,32%, a w porów- 
naniu z najgorszym typoszeregiem T3 zmniejszyła się o 10,85%. 


POLIOPTYMALIZACJA TYPOSZEREGU PRZEKROJÓW 


Przeprowadzona analiza optymalizacyjna dała odpowiedź na pytanie, 
który z typoszeregów minimalizuje wskaźnik zużycia materiału. Nie rozstrzyg- 
nęła natomiast poruszonego na wstępie problemu liczebności typoszeregu. 
Kryteria zestawione w funkcji (19) są w stosunku do siebie konfliktowe. 
Poprawa (zmfiejszenie wartości) jednego kryterium powoduje pogorszenie 
(zwiększenie wartości) drugiego i odwrotnie. Rozstrzygnięcie tego konfliktu 
można uzyskać przeprowadzając analizę polioptymalizacyjną [22, 23]. Polega 
ona na uporządkowaniu wektora funkcji celu w taki sposób, aby poszczególne 
jego składowe nie mogły być polepszone bez jednoczesnego pogorszenia 
pozostałych. Uporządkowanie to prowadzi do wyznaczania rozwiązań nie- 
zdominowanych, czyli Pareto-optymalnych. Mając do dyspozycji zbiór kom- 
promisów można, po uwzględnieniu dodatkowych wymagań i informacji, 
wybrać typoszereg preferowany. 

Wektor x, jest rozwiązaniem zadania w sensie Pareto (należy do zbioru 
rozwiązań kompromisowych ZX), jeżeli odpowiadający mu wektor funkcji 
celu f, (x,) jest wektorem najmniejszym w sensie częściowego uporządko- 
wania oraz jeżeli f, (x,) należy do łamanej wyznaczającej brzeg obszaru celu 
Yi łączącej punkty narożne określone przez minima składowych wektora 
funkcji celu [24]. 

Sposób wyznaczania dyskretnego zbioru kompromisów uzależniony jest od 
charakteru zmiennych decyzyjnych. Jeżeli choć jedna z nich powoduje wielo- 
modalną zmienność funkcji celu, to uzyskanie pełnego zbioru rozwiązań 
kompromisowych ZX i zbioru ocen rozwiązań kompromisowych ZY możliwe 
jest po przeanalizowaniu całego zbioru ocen Y[24]. W rozpatrywanym 
przypadku zmienna T' ma taki charakter, należy więc zaprojektować i ocenić 24 
kratownice. 

Wyniki analizy funkcji y (x,) i ó„,,(%,) pokazano na ryc. 10. W celu lepszego 
przedstawienia zmienności funkcji narysowano ich warstwice. Są one wy- 
znaczone na podstawie dyskretnych wartości funkcji w punktach x,. Trzecia 
z rozpatrywanych funkcji J (x,) zależy tylko od zmiennej k. Wartości J (x,) 
spełniają ograniczenie (22). 

Wyniki dwukryterialnej optymalizacji funkcji (19) pokazano na ryc. 11. 
W obszarze dopuszczalnym X linią łamaną połączono efektywne wartości 
zmiennych decyzyjnych należące do zbioru rozwiązań kompromisowych 
ZX=(xg, Xg, Xj3, X22]. Zbiór ten wyznaczony jest na podstawie analizy 
zbioru ocen Y. 
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Ryc. 10. Warstwice wskaźnika zużycia materiału (a) i przemieszczenia środkowego węzła 
kratownicy (b) 


Rozwiązanie zadania polioptymalizacji (19) wykazuje, że typoszereg dający 
najmniejszy wskaźnik zużycia materiału zmienia się w zależności od liczby 
przekrojów J. Typoszereg T6 najlepszy w przypadku J=25 ma najwięcej 
przekrojów A,<11 cm”, wybieranych w 70—80% przypadków. Typoszereg T1 
najlepszy dla J=13 i J=7 ma tak dobrane przekroje (4,, ,=4 Aj), że pozwala 
na projektowanie kratownic o najbardziej wyrównanych naprężeniach. Typo- 
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Ryc. 11. Zbiór rozwiązań kompromisowych (a) i ocen kompromisowych (b) zadania (19) 


7 Szczecińskie Rocz. Nauk. VI, z. 1 
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szeregi T1 dla J=7i T4 dla J=4 mają najwięcej przekrojów o powierzchni 
A;<11 cm”, co przy danych parametrach L i P ma decydujące znaczenie 
w minimalizacji wskaźnika zużycia materiału. 

Wyniki polioptymalizacji wektorowej funkcji (20) pokazano na ryc. 12. 
Składowymi wektora (29) są wskaźniki zużycia materiału y(T) oraz przemiesz- 
czenie środkowego węzła kratownicy 0,,,(T). Funkcje te tworzą szeroki kąt 
kooperacji a (ryc. 12). Wektor rozwiązań kompromisowych składa się z 19 
kratownic. Tylko 5 ocen rozwiązań nie wchodzi do zbioru ocen kom- 
promisowych. Z punktu widzenia projektanta rozwiązanie to jest niewłaściwe, 
gdyż utworzony zbiór kompromisów pokrywający 79% obszaru dopuszczal- 
nego utrudnia wybór rozwiązania preferowanego. Mimo iż kryteria y(T) 
i 0„„,(7T) są konfliktowe — minima ich leżą na przeciwległych ograniczeniach 
(ryc. 10), to wartości funkcji są ze sobą ściśle powiązane. Wzrostowi masy 
konstrukcji, a więc jej sztywności, towarzyszy proporcjonalny spadek przemie- 
szczeń węzła [8]: 

ŚqL* 


Ó ń 
aD 


(31) 


gdzie 

q — obciążenie M/m?, 

Ę — współczynnik zależny od sposobu podparcia, 
D — sztywność zgięciowa kratownicy. 


Właściwość ta jest szczególnie widoczna na ryc. 10. Warstwice funkcji y(T) 
i ó„.(T) mają zbliżony przebieg, różnią się zwrotem wektora gradientu. 


ROZWIĄZANIE PREFEROWANE 


Rozwiązanie preferowane określono dla wektora f1, (x,) (19). Przy jego 
wyborze zastosowano metodę funkcji dystansowych oraz metodę kryterium 
globalnego. 

Metoda funkcji dystansowych polega na określeniu odległości punktu 
idealnego f,, od poszczególnych punktów należących do zbioru ocen kom- 
promisowych ZY. Współrzędne punktu idealnego wyznaczone są przez mini- 
malne wartości funkcji y(x,) oraz J(x,) (ryc. 11). 

Funkcje dystansowe mają postać [25] 


1 

F*=t> (e) —v (£,)] + [J (e) —J(Ż,)]*)e, ce (1, 2, 0). (32) 
Dla c=2 wartość FÓ) oznacza odległość geometryczną punktu idealnego od 
ocen rozwiązań kompromisowych. Minimalną odległość F'? =3,962 wskazują- 
cą rozwiązanie preferowane zaznaczono na rycinie 11. Zestawienie wartości 
funkcji dystansowych dla uzyskanego zbioru kompromisów zawiera tab. 3. 
Niezależnie od zastosowanej miary odległości metoda funkcji dystansowych 
wskazuje rozwiązanie preferowane dla p=13. Jest to typoszereg T1 o liczbie 
przekrojów J=7. 
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Ryc. 12. Zbiór rozwiązań kompromisowych (a) i ocen kompromisowych (b) zadania (29) 


Wadą tak określonego rozwiązania preferowanego jest jego zależność od 
przyjętych jednostek funkcji y(x,). Przy odpowiednim doborze jednostek 
funkcji celu metoda funkcji dystansowych może wskazać jako rozwiązanie 
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preferowane każdy z wektorów należących do zbioru ocen rozwiązań kom- 
promisowych ZY. 

W metodzie kryterium globalnego przeanalizowano funkcję kosztu porów- 
nawczego konstrukcji. Uwzględniono koszt materiału oraz wykonania kratow- 
nicy. Dla rozpatrywanych zmiennych decyzyjnych założono, że koszt materiału 
jest proporcjonalny do liczby przekrojów w typoszeregu. Globalna funkcja celu 
ma postać 


—, YO) | J(x7) 
FGtx)=", „ej +t2 (zyj PE (6,7, 13, 22). (33) 


Jako porównawcze przyjęto rozwiązanie x,=[6, 25]. Współczynniki wag 
funkcji globalnej są unormowane (c; +c,=1,0), czyli FG(x,)= 1,0. 


|F6 (Xp) 


4.4 


4.2 Ę 


0.6 


0.1 


0.2 


Ryc. 13. Globalna funkcja celu dla różnych współczynników wag 


Podobnie jak w analizie funkcji dystansowych każde z rozwiązań kom- 
promisowych może być rozwiązaniem preferowanym w zależności od wartości 
współczynników wag. Dla c<0,35 preferowane jest rozwiązanie p=22, dla 
c, e 40,36, 0,69) rozwiązanie p=13, dla c, e 40,69, 0,91) rozwiązanie p=7, 
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dla c, 20,92 rozwiązanie p=6 (tab. 3). Przebieg funkcji FG(x,) dla różnych 
wartości współczynnika c, pokazano na rycinie 13. Projektant posiadając 
informacje o proporcji kosztu materiału do kosztu wykonania w danych 
warunkach, może określić wielkości współczynników wag c, ic, i jednoznacz- 
nie wskazać rozwiązanie preferowane. 

Rozkład sztywności prętów dla kratownicy p=13 (typoszereg T1 o liczeb- 
ności j=7) pokazano na ryc. 1. Przekroje prętów opisano przez podanie dla 
każdego z nich wartości j. 


UWAGI KOŃCOWE 


W dyskretnej optymalizacji konstrukcji prętowych istotne znaczenie ma 
właściwy dobór geometrii przekrojów. Budując odpowiednio typoszereg prę- 
tów o określonej liczebności można uzyskać oszczędności zużycia materiału 
rzędu 10%. 

Poprawnie skonstruowany typoszereg przekrojów o ograniczonej liczebno- 
ści powinien w przypadku kratownic spełniać następujące zalecenia: 

— profile należy wybierać z aktualnego katalogu programu produkcji wyro- 
bów, 

— pole przekroju poprzecznego kolejnych elementów typoszeregu powinno 
wzrastać zgodnie z formułą A,,,=4*4j;, gdzie A e (1,1; 1,8), 

— wzrostowi pola przekroju musi towarzyszyć równomierny wzrost siły 
krytycznej, 

— z określonej wielkości kszałtowania należy wybierać taki, który posiada 
maksymalny promień bezwładności przekroju, 

— liczba przekrojów powinna wynikać z kompromisu pomiędzy stopniem 
wykorzystania naprężeń, a kosztem wykonania konstrukcji. 

Ostatnie zalecenie może być realizowane w dwojaki sposób. Liczba różnych 
przekrojów zastosowanych w konstrukcji jest ograniczana przez zmianę 
liczebności typoszeregu przekrojów lub przez wprowadzenie stref jednakowej 
sztywności prętów. W pierwszym przypadku zmienną decyzyjną jest liczebność 
typoszeregu, w drugim zaś zakres i liczba stref sztywności konstrukcji. 
Strefowanie sztywności kratownic przestrzennych często stosowane w prak- 
tycznych realizacjach [18, 21] oraz w licznych pracach teoretycznych [1, 2, 16, 
26, 27] podjęte zostało w aspekcie problemu polioptymalizacji w pracach 
[28, 29]. 


WITOLD MAREK PACZKOWSKI 


THE BARS TYPE RANGE DISCRETE POLYOPTIMIZATION 
OF SPACE TRUSSES 


Summary 


The discrete polyoptimization problem of the ortogonal space trusses point supported in 
corners is formulated and solved. Two variables are considered — the cross section type range of the 
bars and number of elements putting out of the main type range. The type range numbers of 
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elements was changed from 4 to 25. The bars mass of the structure, middle joint displacement and 
the numbers of the bars type range are taken as the objective functions. 

The compromise solutions set is determinated by the systematic search method. The 
preferated solution is determinated by the global criterion method. 

The aim of the work is describing the limited numbers type range usefull in optimal space 
trusses building. 


BHTOJIBJĄ M. IIAUKOBCKH 


NUCKPETHAA IIOJIMOIITAMKH3ANM4A THMIIOBOTO PAJIA 
CTEPKHEMA IIPOCTPAHCTBEHHOŃ O©EPMBbI 


Pe3rome 


B paóore cjopMyjmpoBaHa m  peliefa NACKpETHaa 3aNada NOJHONTUMHAZANAM 
OPTOTOHAJILHOH NpocTpaHcTBeHHOń (hepMBI Ha  HETRIDEX  TOHEYHKIX OIIOPAX B yTJIaX. 
IIpoaHa.1a3ApDOBAHKI J[Be pELIAKOINHE NepeM€HHBIE THIIOBOTO PAĄA NOMepe1HBIX CedeHHii 
cTepxHeń, Aa TaKXKE H4ACJIO 3JIEMEHTOB, BEIAeDKABAEMBIX H3 OCHOBHOTO THIIOBOTO PANA. 
IIocTpoeHbl 6 THMOBBIX pANOB C ACIOJIb3OBAHHEM KATAJIOTA TOpA1eKATAHHKIX CTAJIbHKIX Tpy6. 
UHCJIO 3JIEMEHTOB THIIOBBIX pANOB A3MEHAJIOCb B IpeneJlax OT 4 qo 25. B ka4ecTBe KpHTEpHeB 
OITAMH3AIIAH NpHHATBI Macca cTepxHeń (bepMBI, repeMelIeHHe  IEHTPAJIŁHOTO  yJJIA 
H 4HCJIEHHOCTb THINOBOTO pAĄa CTEepxHeii. MHOXEcTBO KOMIIDOMACCHBIX peLIeHHH 3ANAA 
NONAONTAMA3ALNAA ONPEJEJJEHO METONOM CHCTEMATH1ECKOTO NACKPETHOTO IIOHCKA AOMNyYCTAMOŃ 
oGnacra. IIpenno4TATEHbHOe peuieAne OHNDPEJEJEHO METOJNOM BAJIOBOTO KPHTEPAA. 

B paó6oTe onpeneJleHbIl NPAHIMAIIBI IOCTPOHKA TANOBOTO PANA NoNepedHbIX CedeHHii 
OTDpaHA1eHHOŃ  AHCJIEHHOCTH,  BO3MOXHBIX  [UIA  HCNHOJIb3OBAHHA HPA  ONTAMAJŁHOM 
(bopMapoBaHHH MpocTpaHCcTBeHHbIX (hepM. 
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O MOŻLIWOŚCI WYKORZYSTANIA 
CZĘSTOTLIWOŚCI RICE'A 
DO OCENY JAKOŚCI REGULACJI 


Instytut Automatyki Przemysłowej Politechniki Szczecińskiej 
ul. Gen. Sikorskiego 37, 70-313 Szczecin 


Streszczenie. W pracy omówiono podstawowe właściwości częstotliwości 
Rice'a w aspekcie zastosowań teoretycznych i praktycznych w układach regulacji 
i sterowania, podlegających oddziaływaniu zakłóceń stochastycznych, dla oceny 
jakości regulacji. Z uwagi na proste możliwości pomiaru częstotliwości Rice'a 
sygnału błędu w trakcie pracy układu regulacji, istnieje możliwość wykorzystania 
tego parametru jako dodatkowej informacji, dla celów adaptacji. W pracy roz- 
ważono częstotliwość Ricea w układzie liniowym jedno- i wielowymiarowym, 
opisanym równaniem stanu. W zakończeniu podano wyniki symulacji komputero- 
wych. 


Słowa kluczowe: charakterystyki statystyczne — częstotliwość Rice'a — Zero 
Crosing Rate — jakość regulacji. 


WSTĘP 


Wraz ze wzrostem wymagań technologicznych, rosną zadania stawiane 
regulatorom przemysłowym. Niezbędne stają się prace nad zastosowaniami 
nowych algorytmów i nowych struktur regulatorów autonomicznych z adapta- 
cją do zmiennych warunków pracy. Niezbędne są także poszukiwania i próby 
tworzenia kryteriów oceny jakości regulacji, uwzględniających losowy charak- 
ter zakłóceń, jakie w rzeczywistości oddziałują na obiekt regulowany i cały 
układ sterowania. 

Rozpatrywanie problemu jakości w kategoriach sygnałów stochastycznych 
jest uzasadnione również względami praktycznymi, aby pomiary wielkości 
charakteryzujących jakość mogły odbywać się w czasie normalnej pracy 
rzeczywistych układów regulacji z rzeczywistymi obiektami. Jednym ze wskażź- 
ników uzupełniających [2] może być częstotliwość Rice'a [1,6-10] sygnału, to 
jest średnia liczba przejść przez poziom zerowy w jednostce czasu, scent- 
rowanego sygnału błędu. Wielkość ta związana jest bezpośrednio ze średnio- 
kwadratową prędkością zmian sygnału błędu i jego wariancją i charak- 
teryzować może „ruchliwość układu” rozumianą jako zdolność do szybkiego 
reagowania na zakłócenia. Praktycznie ważne jest też i to, że może ona być 
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prosto mierzona w trakcie normalnej pracy układu regulacji. Większość 
przemysłowych procesów stochastycznych w systemach sterowania automaty- 
cznego i regulacji spełnia teoretycznie wymagania nieograniczonego różnicz- 
kowania, jednak ze względów praktycznych, z uwagi na ograniczoną rozdziel- 
czość przyrządów pomiarowych, praktycznie można oceniać częstotliwość 
Rice'a sygnału lub jego pochodnej, to jest można mierzyć praktycznie 
w zasadzie tylko średnią liczbę przejść sygnału przez zero lub określony 
poziom oraz mierzyć średnią liczbę maksimów sygnału, przypadających na 
jednostkę czasu. Można wykazać, że znając częstotliwości Rice'a sygnału i jego 
pochodnych, uzyskuje się informację równorzędną z tą, jaką daje funkcja 
autokorelacji tego sygnału R(z) dla niewielkich wartości r [8, 9, 10]. 


CZĘSTOTLIWOŚĆ RICEA W UKŁADZIE REGULACJI 


Rozważymy stacjonarny, liniowy układ zamknięty, jak na ryc. 1. Załóżmy, 
że układ pracuje w reżimie stabilizacji (xy(t)=const) i poddany jest za- 
kłócającemu działaniu sygnału stochastycznego z (t) od strony wejścia obiektu. 
Przyjmijmy, że z (t) jest sygnałem stacjonarnym scentrowanym (E (z(t))=0), 
charakteryzowanym autokorelacją R.(r), a obiekt i regulator są opisane 
transmitancjami, odpowiednio K,(s) i KZ(5). 


Obiekt 
Kgls] 


z(t) 


Regulator 
Kąę (s) 


KE = ZZOZ RÓ ZZ ZA] 


Ryc. 1. Schemat blokowy rozważanego układu 


O jakości regulacji można sądzić w kategoriach sygnałów stochastycznych 
na podstawie dwu podstawowych parametrów [2, 3] 
— wariancji sygnału błędu R„(0), 
— średniokwadratowej szybkości zmian sygnału błędu 


= (3) 


Średniokwadratowa szybkość zmian sygnału błędu, charakteryzuje szybkość 
działania układu regulacji i jest związana z częstotliwością Rice'a sygnału błędu 
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e(t), która ogólnie, dla stacjonarnego, scentrowanego sygnału e(t) o auto- 
korelacji R,(r) jest określona wzorem 
17. R; BEM 
4,=— lim A e” 2R.(0) , (1) 
T +0 /1—R32(2) 
gdzie c jest poziomem, dla którego określa się częstotliwość Ricea. 
Przy czym, dla przejść przez zero (c=0) i dla różniczkowalnych średniokwa- 
dratowo procesów e(t), otrzymuje się w granicy na podstawie (1) zależność 


_1 [-RE0) 
tez PRO) » 
gdzie 

nę _6%R.(0)| _„ffde(y? 

R0-77"|„„zfa) |: s 
Ponieważ 

R,(j=-- | 6,(ajeś do, (4) 
21 


więc po dwukrotnym różniczkowaniu, dla 1=0, będzie 


co 


= . | w*G,(0) do. 6) 
21 


—o 


_62R,(2) 


a1? 


1=0 


Stąd zależności (2) w dziedzinie czasu odpowiada równoważna postać w dzie- 
dzinie częstotliwości 


[ 0*G,(0)do 
o |= 


e 


pe = , (6) 
| G,(o)do 


gdzie G,(w) jest gęstością widma mocy sygnału błędu: 
G,(0)=|K (jo)|? G. (0). (7) 


Zatem zależność (6) może być przedstawiona w wygodniejszej postaci 


1 [ w?|K (jo)|?G, (0) dv 
0 


1ę= (8) 


nż © 
[ |Kdo)?G, (0) do 
[o 
gdzie, zgodnie z ryc. 1 

K,(iw) 
1+K, (jo) Kz(jo) 


K (j0)= ©) 
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Rozważając problem w dziedzinie czasu można napisać 
t © 
e(t)= | k(t—o)z(a)do= | k(o)z(t—o)do, (10) 
- 0 [U 


gdzie k(a) jest impulsową funkcją przejścia układu z ryc. 1. Różniczkując (10) 
będzie 


e(t)= | k'(t—a)z(a)da +k(0)z(t). (11) 


Stąd można napisać 
R. (t, t)=Ete(t)e'(t,)).. (12) 


Stosując podstawienie r=t, —t i odpowiednio przekształcając uzyskuje się 


R: ()= | | k(o)k(ty R.c-0+3)dody+k(0) [ K)R.(r+7)dy+ 
00 [o] 
+k(0) | k(o) R.(6-2)do+k?O)R„(a). (13) 
[e] 


Uwzględniając (13) we wzorze (2) dla r=0 można napisać 


0 


1 | [OKR (-odody+k(0 [EGJR,0)óy+KO [KOR Gab +ŁLOR,() 
a 


q 0 © 

| | Kok()RZÓ o) dydo (14) 
Jest to postać czasowa wzoru (8). Ze względu na swe skomplikowanie, 
w porównaniu z (8), wzór ten może okazać się nieefektywny. Ponadto należy 
zauważyć, iż o ile realne procesy stochastyczne są gładkie z natury i dla nich 
wzór (14) jest w pełni słuszny, o tyle w przypadkach abstrakcyjnych modeli 
procesów w postaci procesów Markowa, z uwagi na ich niegładkość, nie może 
on być stosowany. 

Okazuje się jednak, iż w takich przypadkach efektywna jest zależność (8), 
z dziedziny częstotliwości, natomiast w dziedzinie czasu należałoby w takich 
przypadkach stosować teorię dystrybucji. 

Można też zauważyć, iż w dziedzinie czasowej z formalnego, matematycz- 
nego punktu widzenia, obowiązują ostrzejsze założenia niż w dziedzinie 
częstotliwości; stąd dla tego samego problemu wzór (14) może się okazać 
nieefektywnym, w przeciwieństwie do zależności (8). Dla przykładu rozważymy 
przypadek, gdy na element inercyjny o stałej T trafia zakłócenie z (t) w postaci 
sumy N sinusoid. Wtedy autokorelacja 


N 
R.(r)=2 4 j cosQ,T. 
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Odpowiada temu gęstość widmowa 


G,(0)=nr y Aj 6(Q,—0). 


j=1 
Stosując zależność (8) uzyskujemy od razu wynik 
N Aj Q 
1 > 1+Q7T2 T2 
CE 26 


2 Traf? t+Q; T* 
Identyczny wynik uzyskuje się (R,(r) gładka w zerze) za pomocą wzoru (14), 
lecz po szeregu przekształceniach. Gdy natomiast przyjąć, dla tego samego 
obiektu, zakłócenie z (t) w postaci sygnału telegraficznego, dla którego auto- 
korelacja jest określona zależnością 
R.(0)=R.(0)e'*"! 


i jest niegładka w zerze, to wzór (14) okazuje się nieefektywny, za pomocą (8) 
zaś, dla 


R(0) a 
G. = z 
(0) n a +w? 
uzyskuje się prosto 
1a 
lż=—3 T . 


W tab. 1 przedstawiono obliczone wartości 42 dla kilku wybranych typowych 
przypadków sygnału wejściowego przyjmując, że charakterystyka całego ukła- 
du może być reprezentowana transmitancją m-inercyjną o jednakowej stałej 
czasowej 

K(s)=k(1+sT) ". (15) 


Dla dowolnego, nieanalitycznego sygnału zakłócającego, którego gęstość 
widma mocy jest znana w postaci wykresu G,(w), np. na podstawie eks- 
perymentu i może być przedstawiona jako suma gęstości procesów pas- 
mowych, jak to przykładowo pokazano na ryc. 2, można na podstawie wzoru 
(8) obliczyć w przybliżeniu AŻ. Przyjmując 0; — 0; _,=4, A=const, G,(0;)=a;, 


dla i=1,..., N, uzyskuje się zależność 
i4 
y wdw 
1 = (Lhor TYP 
A = m 4 : | (16) 
y do 
= JS (FET 
(i- 1) 4 


co uwzględniono również w tab. 1. 
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Gz (w) 
j |] 


Wi 
OB; _ą 


Ryc. 2. Dyskretyzacja charakterystyki G,(w) dla obliczeń przybliżonych 42 


W tab. 2 zestawiono kilka wyników obliczeń AŻ dla układu z ryc. 1, 
w przypadku obiektu m-inercyjnego o stałej T oraz regulatora PID o trans- 
mitancji Kp(s)=Kz;(1+STą+1/sT,). We wzorach przyjęto k=K,K,. 

Tak więc mierząc lub dokonując analizy A, w układzie regulacji można 
uzyskać informacje o zmianie sytuacji w układzie z ryc. 1. W przypadku 
pomiaru, informacje te mogą być.uzyskiwane i wykorzystywane na bieżąco dla 
celów adaptacji lub dla innych celów. Można zaproponować, podobnie jak 
w [2, 3], wskaźnik dający obraz zwiększenia „ruchliwości układu” w postaci 


AŻ z reg. 


Aż bez reg. s 


dą 
Wskaźnik ten, uwzględniając w (17) wzory (6) i (9), dla sytuacji jak na ryc. 1, 
można przedstawić w postaci zależności 


0 
2 


K;, (jw) 
JG 2 0 
Jk. a |= |rztokio| 94a 
0 
K, (jo) ; 
|= o (0) "6, (w)do Jo*KoboG,w)da 
ć (18) 
lub odpowiednio w dziedzinie czasowej 
s [(G7) | rs 
z reg. 
? d 
_ R.(0) bez reg (19) 


GR Ol2Gi "zWZD(ANZE = == 6 
R.(0 z reg. (DĄ SE 
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Zatem wskaźnik q, można interpretować jako multiplikatywny wskaźnik 
„wielokryterialny”, uwzględniający jednocześnie redukcję wariancji (pierwszy 
człon prawej strony równań (18) i (19)) oraz zwiększenie średniokwadratowej 
szybkości sygnału błędu (drugi człon prawej strony), dzięki zastosowaniu 
regulatora. Można łatwo obliczyć (18) dla przypadku, gdy zakłócenie działają- 
ce na układ od strony wejścia obiektu jest sumą N sinusoid. Wtedy będzie 
N N 202 o 2 
Y ałlkoaj? © roewiwkak=, 
k=1 za 1+Ko (60) Kz(9,)| 
__ AEJKoGO.J? | 


N 
EE M NERSĄOGZ POP OMIŻ A? QŻ K. GQ 2 
 FK,GQJKAGOE 2, oi) 


W szczególnym przypadku, gdy N=l, q,=l1, niezależnie od postaci trans- 
mitancji obiektu i regulatora. Oznacza to, iż w przypadku szczególnego 
zakłócenia w postaci pojedynczej sinusoidy, wskaźnik q, nie jest wrażliwy na 
zmiany właściwości dynamicznych układu regulacji. 

Opierając się na spostrzeżeniu, iż mierząc 4, oraz 4, gdzie 4,, jest 
częstotliwością Rice'a dla poziomu cz40, można pośrednio mierzyć i szacować 
wariancję gaussowskiego procesu e(t), jeżeli spełnia on założenie o normalno- 
ści. Wtedy, na podstawie (1), jest 


4,= 


Jz 


k 


II 


cz 


R,(0)= (20) 


y= 
Żln 
A 


Zatem, można mierzyć w ten sposób wariancję błędu e(t), a także można 
wykorzystać drugi wskaźnik, który mógłby być określony na podstawie 
pomiaru częstotliwości Rice'a, dający obraz zmniejszenia wariancji dzięki 
zastosowaniu regulatora lub też dzięki odpowiedniemu doborowi jego para- 
metrów. Byłby to iloraz wariancji błędu w układzie z i bez regulatora. 


n 1 bez re 
— R.(0) z rcg. <zj A Ę | (21) 
de" R_(0) bez reg. (>) 
In "a z reg. 


Stosując wskaźniki q;, qz można porównywać różne algorytmy w obecności 
zakłóceń stochastycznych, można również ocenić i porównać „dzielność” 
regulatorów pracujących w obecności zakłóceń losowych [3] z określonym 
obiektem. Mierząc A, 4,, można wykrywać zmiany parametrów obiektu 
w czasie jego pracy lub też zmiany parametrów i właściwości zakłóceń. 
Przyjmując ponadto, że c jest dopuszczalnym poziomem, określającym pas 
tolerancji dla sygnału błędu, można analizować częstotliwość przejścia tego 
poziomu lub średni czas występowania takiej sytuacji w układzie regulacji. 

Przykładowo, dla obiektu m-inercyjnego przy m=3 i regulatora PID (tab. 
2) przy zakłóceniu o autokorelacji R,(1)=e "Fl dla przyjętych wartości 


g* 
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parametrów «T=1, Tą=0, 1T, T,=T, uzyskuje się dla różnych wzmocnień 
w pętli następujące wartości wskaźnika (17): 


bez reg. q,=1 
k=l 41 =5 
k=9 q,=25. 


CZĘSTOTLIWOŚĆ RICEA W LINIOWYM UKŁADZIE 
WIELOWYMIAROWYM 


Rozważymy układ stacjonarny liniowy n-tego rzędu o p wejściach i m wyj- 
ściach, opisany równaniem stanu i równaniem wyjścia: 


> » 
gdzie macierze A, B, C i wektory x, u, y są odpowiednio o wymiarach 
A:nXn, 
B:n xp, 
C:mxn, 
x:nxl, (23) 
u:p xl, 
y:mxl. 
Dla i-tego wyjścia, gdzie i= 1, 2, ..., m, kwadrat częstotliwości Rice'a A, będzie 
„2 
A (24) 
gdzie uwzględniając (22), ogólnie 
k = Z. + Bu) ? 


Tworzymy macierz H o wymiarze mxm będącą iloczynem zewnętrznym 
wektora p 


H=jj"=(hyjlnxu* (26) 
Analogicznie tworzymy macierz G o wymiarze mx m dla wektora y 
G=yy"=lgylmxm* (27) 


Skąd, uwzględniając (24), określimy i-ty kwadrat częstotliwości Rice'a w po- 
staci ilorazu wartości oczekiwanych, zgodnie z poniższą zależnością 
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n? Etg) 
gdzie h,,, g,, są elementami diagonalnymi macierzy, odpowiednio H oraz G. Na 
podstawie (26) i (27), uwzględniając (22) można napisać 


H=C(Ax + Bu)[C(Ax+ Bu)]7, 
H=C(Ax+ Bu) (Ax + Bu)" CT, (29) 
G=CxxTCT. (30) 


Elementy h,, oraz g,;, po żmudnych przekształceniach, określone są zależnoś- 
ciami 


hu 2 X%% 2 Ak Ci 2 a,,C, + 
cf 5 5 UĄ U, > by Ci, » bc, +2 » XK U > Ak Cię >, Być ir > (31) 


k=1l=1 k=1 l=1 
EBY Ci X; s. (32) 
r=1l 


Obliczając na podstawie (31) wartość oczekiwaną E £h,,), uzyskuje się wyraże- 
nie ogólne 


Eth,) = — > > Rz, xi (0) l > App Cir” > a, c] cj 
k=1 I= r= r= 


1 


Pp p n n 
PZ Z Rad! > bik Cie” ą byc] + 
k=1 I=1 r=l1 r= 
n p 
+2 >. > R uzl OJ9: Azk Cir > b4ćj,)- (33) 
k=1 I=1 


Tak więc licznik wyrażenia (28), określającego kwadrat częstotliwości Rice'a 
dla i-tego wyjścia, jest związany z korelacjami w zerze poszczególnych 
zmiennych stanu, korelacjami sterowań w zerze oraz z korelacjami zmiennych 
stanu i sterowań w zerze, zaś wyrażenia w nawiasach kwadratowych są 
współczynnikami liczbowymi. Dla układów z jednym wejściem i jednym 
wyjściem (SISO) m=1, p=1, wzory się znacznie upraszczają 


Elhd=Y YR..OLY asc, X 4,6] + 
r=l r= 


k=1 I=1 


[X bc]? 


+22 Rz( OCZ ane, be. (34) 
1 
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Mianownik wyrażenia (28) będzie w tym przypadku określony zależnością 
E(gul=L > Rx, (0) CC, . (35) 
k=1 1I=1 


Jeżeli rozważany układ opisany jest transmitancją 
by+b,S+...+b,s" 
ay+a, S+...+a,_1+57" 


to stosując postać kanoniczną fazową można ten układ opisać równaniem 
stanu i równaniem wyjścia w postaci 


K(s)= (36) 


x, (0) 6: 0£ 062 0 x, (1) 0 
*.() = 0 0 1 sa 0 *z() 4 0 U(), 87) 
X, (t) —09 —044 —0 ... —A,-Y x„(t) 1 
y(t)=[bob; -.. b, 0...0] [ x; (t) 
AA (38) 
x,(0 


Rozpatrując ogólnie SISO, z różniczkowaniem v-krotnym w (36) i stosując 
oznaczenia dla elementów macierzy c 


6 =B; 
Cz =by 
cą =b, 
q=v+1 


można po przekształceniach napisać 
E (911) =R;, xi (0) d+R,., x2 (0) lą G+R;, x3 (0)c, C3 RR 94 xq () C1 cat 
z: Rz, x1 (0) cą Cą + Rz, x2 (0) cą ca FRZ, x3 (0) cą c3 +... Ra Xxq (0) c (0)7 (39) 
E LO, =(1 KA Baa (0)+ Ca Ry (0)+ C3 AAC „Tt cą R...(0)] k 
A [GR zt „(0)+c2 Raga (0) + cz Ra „(0)+...+ ć, Bo: ŻR* (0)]. (40) 


Ostatecznie na podstawie powyższych zależności można napisać ogólnie wzór, 
określający kwadrat częstotliwości Rice'a dla układu SISO w postaci kanonicz- 
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nej, n-inercyjnego z v-krotnym różniczkowaniem 


q q 

1 2 C; > c,R;,,, szl) 

42 =— iF1_ 351 
nr ą ą ” 
Ż © 2 GR... (0) 
i=1 j=l 

Wzór (41) można przekształcić do nieco innej postaci, aby wyróżnić człony 

z autokorelacją w zerze oraz człony z korelacją wzajemną w zerze współrzęd- 

nych wektora X 


(41) 


ą 4-1 4 
1 pa cj Ri, x, .,(00 +2 > c; DRZCZ 0) 
a= RansdaB SR 07 (42) 
> c? Ryl Xi 0)+2 paz PAR R, «,(0) 
i=1 j=i+1 


Zależności (41) lub (42) ilustrują Ads przedstawienia kwadratu częstot- 
liwości Rice'a w postaci ilorazu kombinacji liniowych autokorelacji i korelacji 
wzajemnych zmiennych stanu dla r=0, aż do q+1 zmiennej włącznie. Tak 
więc rząd różniczkowania w liczniku transmitancji (36) decyduje o liczbie 
kombinacji. Ponieważ q=v+1, więc jeśli rząd różniczkowania jest v, to 
korelacje występują, aż do v+2 zmiennej stanu. 

W praktyce, z uwagi na kłopoty z różniczkowaniem sygnałów zakłóconych 
szumami pomiarowymi, nie stosuje się z reguły wyższego różniczkowania niż 
v=2. Najczęściej (np. regulator PID) stosowane jest jednokrotne różnicz- 
kowanie (v=1). Stąd w takich przypadkach wzory, określające kwadrat 
częstotliwości Rice'a, będą stosunkowo proste. Jeżeli transmitancja (36) odnosi 
się do układu z ryc. 1 to wyraża ona stosunek transformat błędu regulacji 
w układzie do wejściowego zakłócenia. 


SYMULACJE 


Dla ilustracji rozważań przeprowadzono symulacje na komputerze IBM 
XT z wykorzystaniem programu „SYMUL”. Program ten służy do roz- 
wiązywania w dziedzinie czasu układów równań różniczkowych pierwszego 
rzędu, połączonych układami równań algebraicznych i napisany jest w języku 
FORTRAN 77. Podprogram dotyczący omawianych symulacji pokazano 
w dodatku 1. Do symulacji wybrano dwie struktury układu regulacji z obiek- 
tem m-inercyjnym i regulatorem PID. W obu przypadkach na układ zamknięty 
działa zakłócenie z(t) o charakterze pseudolosowym, przy czym w przypadku 
a) zakłócenie to działa na wejście obiektu symulując, przykładowo, wahania 
napięcia zasilania, zaś w przypadku b) zakłócenie z(t) symuluje szumy 
pomiarowe i działa na wyjściu obiektu. Te dwie struktury można też 
interpretować jako ilustrację eksperymentu czynnego (a) i biernego (b) w przy- 
padku identyfikacji zmian parametrów obiektu, regulatora czy też zmian 
warunków pracy tego układu. 
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Dla celów omawianych symulacji przyjęto m=4, zaś nastawy regulatora 
PID dobrano wg kryterium optymalnego modułu, na podstawie przytoczonej 
w [5] tablicy dla m=4 i stałej czasowej obiektu T=1 s. Parametry badanego 
układu w postaci rzędu mianownika transmitancji K,(s) lub K,(s) oraz 
współczynników wielomianów licznika i mianownika tych transmitancji wpro- 
wadzono do podprogramu symulacji po doprowadzeniu transmitancji do 
postaci (36). Symulacji dokonano dla stałej czasowej obiektu T=0,8 s, 1 s 


AWIW 


Ryc. 3. Przykładowe przebiegi sygnałów z(t) oraz y(t) dla struktury K,(s), przy różnych 
wartościach stałej czasowej obiektu, dla m=4 
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0 1000 s 


Ryc. 4. Liczba przejść przez zero sygnału z(t) w funkcji czasu, dla symulacji z ryciny 3 


i 1,5 s, nie zmieniając nastaw regulatora PID. Pozwoliło to na zilustrowanie 
wpływu zmian parametrów na częstotliwość Rice'a na wyjściu układu. 

W wyniku symulacji uzyskano wartości wariancji i kowariancji zmiennych 
stanu układu, liczbę przejść przez zero dla sygnału wejściowego i wyjściowego 
oraz wykresy przebiegów w czasie 200 s. 

W tab. 3 przedstawiono niektóre wyniki symulacji i przytoczono dane 
z wydruków, które pozwoliły sprawdzić słuszność wzorów (41) i (42). Na ryc. 
3 pokazano przykładowe przebiegi z(t) i y(t), uzyskane w wyniku symulacji dla 
struktury K,(s) i A,=0,47 1/s, przy różnych wartościach stałych czasowych 
obiektu: T=0,8 s przebieg b, T= 1,0 s przebieg ci T=1,5 s przebieg d. Na ryc. 4 
pokazano liczbę przejść przez zero sygnału z(t), jako funkcję czasu symulacji. 
Na ryc. 5—10 pokazano wyniki symulacji przy zmianie A, w zakresie od 0,02 
1/s do 2,5 1/s. 


UWAGI KOŃCOWE 


Przeprowadzone dla przypadku układu liniowego i stacjonarnego roz- 
ważania teoretyczne pozwalają stwierdzić, że częstotliwość Rice'a może być 
dobrym, uzupełniającym parametrem, charakteryzującym pod kątem jakoś- 
ciowym dynamikę systemu automatyki w kategoriach procesów stochastycz- 
nych. Z uwagi na prostotę pomiaru tego parametru w trakcie działania 
systemu automatyki, możliwe jest praktyczne wykorzystanie tego pomiaru do 
bieżącej oceny jakości regulacji, na podstawie zaproponowanych w pracy 
wskaźników, traktowanych jako uzupełniające do innych kryteriów oceny 
jakości jak na przykład wariancji sygnału błędu w układzie. Dzięki wspo- 
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mnianej prostocie pomiaru częstotliwości Rice'a na bieżąco, wynik tego 
pomiaru może być wykorzystany w praktyce dla celów adaptacji regulatora. 
Może to być zrealizowane na drodze odpowiedniego uzależnienia parametrów 
algorytmu sterowania od wartości częstotliwości Rice'a błędu regulacji i ewen- 
tualnie jego pochodnej. Pozwala to na konstrukcję prostych algorytmów 
adaptacyjnych. 

Wyniki symulacji dobrze ilustrują wrażliwość częstotliwości Rice'a na 
zmiany warunków pracy układu regulacji, czy też na zmiany parametrów 
dynamicznych tego układu. Jednocześnie wyniki symulacji wskazują na to, że 
ocena jakości regulacji przeprowadzona tylko na podstawie pomiaru czy też 
analizy A,, może okazać się niewystarczająca i winna być uzupełniona 
pomiarem czy też analizą innych kryteriów. Symulacje potwierdzają pogląd, iż 
wzrost ruchliwości układu jest w procesie regulacji korzystny z punktu 
widzenia jakości regulacji. Wynika to z faktu, iż w paśmie regulacji zmniej- 
szeniu wariancji błędu regulacji towarzyszy wzrost częstotliwości Rice'a tego 
błędu. Ograniczenia zastosowania mogą wynikać z tego, że częstotliwość 
Rice'a, podobnie zresztą jak i wariancja, jest tylko parametrem liczbowym, 
choć jak wykazano jest ona w sposób bardzo złożony zależna od parametrów 
układu i cech statystycznych zakłóceń działających na układ, charakteryzując 
średniokwadratową szybkość zmian sygnału błędu. Przyjmując błąd i jego 
pochodne jako stan systemu dynamicznego, można analizować zatem średnio- 
kwadratową szybkość zmian stanu układu. Jednak o ile minimalizacja warian- 
cji błędu regulacji, gwarantuje zachowanie odpowiedniej jakości sterowania, 
o tyle maksymalizacja średniokwadratowej szybkości zmian stanu układu, czy 
też maksymalizacja częstotliwości Rice'a A, tego nie może zagwarantować. 
Stąd też wynika charakter uzupełniający tego wskaźnika. 
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Ryc. 8. Porównanie charakterystyk o? 
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=4 


1,0 s przy zmianie 7 obiektu. Struktura K,(s), m 


Ryc. 9. Charakterystyki 2, =o(2,) uzyskane w wyniku symulacji dla obiektu z regulatorem PID, nastrojonym wg kryterium optymalnego modułu, dla 
T: 
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op AIS ISoupo $ EĄAISAIS1ĄEIEYJ SGT=J ISQT=L S80=/L wolsoem tfepermodpo '€ I 7 *T INAISAJS]ĄEIEUYD :SQT=J €£=W ejp siujewAjdo wAuofoljseu 
Qld walojejndal z (s)"y sz1mĄnIs o OdaJ3IUĄWEZ NpEłĄN EJp ZEJO u YoAuzo1 Kzid niysiqo O3awes ejp (7)0©=77 ĄAKISAISIĄEIEYD SIUŁUMOJOJ '0T 249 


Wj | 


; l-5]3Y 
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$STORAGE: 2 

$FLOATCALLS 

EG ZR JEŻA EZ ZZZZZZZZ ZZ K 
C** EQUA 


SUBROUTINE EQUA(T,X,P) 

REAL*4 T,X(1).P(1) 

REAL*4 CON(20).TB,PAR(1000) ,XX(500) „TB1 

REAL*8 DT 

INTEGER*2 ISTR(10),N,N1,N2,1.N3,NU,J,K 

COMMON CON,TB,DT,XX,PAR,ISTR 
(8 
C** Podprogram umożliwia symulacje skalarnego układu dynamicznego 
C** równań stanu w postaci kanonicznej sterowalnej n-tego rzędu 


«* Parametry * 


PAR(1) - stopień mianownika transmitancji 

PAR (2) -— współczynnik a0 mianownika 

PAR(n+1) — współczynnik an-1 mianownika 

PAR (n+2) — współczynnik b0 licznika 

PAR(2n+1) — współczynnik bn-1 licznika 

PAR(2n+2) — parametr sterujący rodzajem sygn. wejściowego 
PAR(2n+3) — wartość początkowa dla programu UNIF "(.12345678901) 


PAR(2n+4) — parametr sterujący odstępem czasu losowania 


* Zmienne całkowane * 


X(1) = 

s — zmienne stanu układu 
X(n) 2 
X(n+1) _ 

: — wartości średnie zmiennych stanu 
X(2n) = 
X(2n+1) — wariancja x1(t)*x1(t) 
X(2n+2) — kowariancja x1(t)*x2(t) 
X(3n) — kowariancja x1(t)*xn(t) 
X(3n+1) — wariancja x2(t)*x2(t) 
X(3n+2) — kowariancja x2(t)*x3(t) 
X(4n-1) — kowariancja x2(t)*xn(t) 
X(4n) — wariancja x3(t)*x3(t) 


NGNDQRGQAGRQQRAQGGGAQRQGNNQGRNQGGNQAQAQGRQQQQQ 


X(4n+1) — kowariancja x3(t)*x4(t) 


(83 — 
[© wzi OŚ A O O A A A w R A A R 
C X(2n+n*(n+1)/2) — wariancja zmiennej stanu xn(t)*xn(t) 
c == NOZNA ZA 
C X(2n+n*(n+1)/2+1) — wartość średnia wyjścia y(t) 
C X(2n+n*(n+1)/2+2) — wartość średnia wejścia u(t) 
c z A RR RZA 
C X(2n+n*(n+1)/2+3) — wariancja y(t) 
C X(2n+n*(n+1)/2+4) — wariancja u(t) 
(6: pam nawao 
E 
(e: EZ O Z YO CYC 
€ « Zmienne dodatkowe * 
(e; 
(e: Z Z WWR RW 
C XX(CON(1)+1) -— wyjście układu y(t) 
C XX(CON(1)+2) — wejście układu u(t) 
e A ZRZEC ZE 
C XX(CON(1)+3) -— liczba przejść dla y(t) 
C XX(CON(1)+4) -— liczba przejść dla u(t) 
EG A A PO RZY 
C XX(CON(1)+5) — zmienna pomocnicza dla programu UNIF 
(© p OP KAPENA 
(e: 
[e i Uwagi * 
[> 
c 
Ę CON(1) = N + N + N(N+1)/2 + 4 - liczba zmiennych całtkowanych 
(e: CON(5) = 4+ 1 — liczba zmiennych dodatkowych 
(e: CON(7) =N+N + - liczba parametrów 
c 
N=INT(PAR(1) ) 
N1=N-1 
N2=N+N+2 
NU=CON (1) +2 
DO 10 I=1,N1 
10 P(I)=X(I+1) 
P(N)=0. 
DO 20 I=1,N 
20 P(N)=P(N)-PAR(I+1)*X(I) 
P(N)=P(N)+XX (NU) 
(© iR EEE 
C** wartości średnie wsp. stanu 
€ z R Z 
N1=N+1 
N2=N+N 


DO 30 I=N1,N2 
P(I)=X(I-N)/CON(3) 
30 CONTINUE 


(8: 
C** wariancje wsp. stanu 
CC NN NN ZEE 
K=N2 
DO 40 I=1,N 
DO 40 J=I,N 
K=K+1 


40 P(K)=X(I)*X(J)/CON(3) 
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C** wartości średnie y(t) i u(t) 
C SSL ccc cw 
N1=K+1 
N2=N1+1 
N3=CON (1) 
DO 50 I=N1,N2 
N3=N3+1 , 
50 P(I1)=XX(N3)/CON(3) 


N1=N2+1 
N2=N1+1 
N3=CON (1) 
DO 60 I=N1,N2 
N3=N3+1 
60 P(I1)=XX(N3) *XX(N3)/CON (3) 

RETURN 
END 
SUBROUTINE EXAM 
REAL*4 CON(20),TB,PAR(1000) ,XX(500),SU,SY,R 
REAL*8 DT,Z i 
INTEGER*2 ISTR(10),N1,N,N2,NX,NS,N4,N3,NY 
COMMON CON,TB,DT,XX,PAR,ISTR 
N=INT(PAR(1)) 
N1=CON(1)+1 
N2=N1+1 
N3=N+N+2 
N4=N3+2 
NX=CON (1) 
NY=N* (N+1)/2+N+N+1 
IF(TB.EQ.0.) THEN 
SY=0. 
SU=1. 
XX(NX+3)=0. 
XX(NX+4)=0. 
Z=PAR(N+N+3) 
NS=-1 
ENDIF 

c — 

C** wyjście układu y(t) 


XX(N1)=0. 
DO 10 I=1,N 

10 XX(N1)=XX(N1)+PAR(I+N+1)*XX(I) 
SY= (SY-XX(NY) ) * (XX(N1)-XX(NY) ) 
IF(SY.LT.0.) XX(NX+3)=XX(NX+3)+1 


SY=XX(N1) 
[0: —u—= muze nn zn NW RRCD CZA 
C** wejście sterujące u(t) 
€ nnn O A A ZO OO 
N1=CON (1) +5 
NS=NS+1 
R=NS,/PAR(N4) 


IF(R-INT(R).EQ.0) THEN 
CALL UNIF(XX(N2),1,Z) 
XX (N2)=XX(N2)-.5 

ENDIF 
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60 IF(PAR(N3) ) 30,40, 20 

20 XX (N2)=PAR(N3) 
GO TO 50 

30 IF(ABS(XX(N2)).GT.0.1) THEN 
XX(N2)=SIGN(PAR(N3) ,XX(N2) ) 
ELSE 
XX(N2)=0., 
ENDIF 

40 SU=SU*XX (N2) 
IF(SU.LT.O. ) XX(NX+4)=XX (NX+4)+1 
SU=XX (N2) 

50 RETURN 
END 
SUBROUTINE UNIF (U, LNR,T) 

c PODPROGRAM GENERUJE LNR CIAGÓW LICZE LOSOWYCH (0,1) 

REAL*4 U(LNR) „CON 
REAL*8 T 
INTEGER*2 IC,I 
CON=3125. 
DO 19 I=1,LNR 
T=CON*T 
IC=DINT(T) 
T=T-IC 
U(I)=T 

10 CONTINUE 
RETURN 
END 


STANISŁAW SKOCZOWSKI 


ABOUT POSIBILITY OF MAKING USE 
OF THE ZERO CROSING RATE 
TO ASSESSEMENT THE QUALITY OF REGULATION 


Summary 


The elementary properties of the Zero Crosing Rate (ZCR) frequency in according with the 
theoretical and practical applications in error-actuated control systems, subordinating of interac- 
tion of the stochastic perturbations, to assessment the quality of regulation has been shown in the 
paper. 

This parameter can be use as an additional information for general adaptation, because the 
measurement of the ZCR of the error signal is very simple during the work of the regulation 
system. 

In this paper has been investigated the ZCR linear system of one or multi-dimensions which 
has been described by the equation of state. 

These results can be use in the computer symulation. 


CTAHKCJIAB CKO4OBCKM 


O BO3MOXHOCTH HCIIOJIB3OBAHHA UACTOT5BI PHNA 
JUIA OIIEHKH KAUECTBA PETYJIMPOBAHHA 


Pe3oMe 


B pa6ore oÓcyK1eHkI TJIaBHkIe CBOŃCTBA WACTOTHI PHINA B ACIIEKTE IMPAKTHŁECKOTO 
HpAMeHeHHA B CHCTEMAX ETyJIApOBAHHA H yIIpaBJIeHHA, IOĄBEPTHYTKIX  BO3NEACTBHKO 
CTOXACTA4ECKHX  IIOMEX, JUIA OLNEHKA KA4ECTBA DpEryJlApoBaHHA. DJIaronapa HpocTkIM 
BO3MOXHOCTAM H3MEpeHHA 4ACTOTRI PHINA CAHTHAJIA OLIAÓKA B XONe paÓOTKI CHCTEMKI 
peryJldpoBaHHA,  CYyNIECTBYET  BO3MOXKHOCTŁ  HCIOJIE3OBAHHA  3TOTO  IapaMeTpa,  Kak 
AONOJHATEJIEHOH HHDOPMALHA, IA neleńi aqauTaqHH. B pa6oTe paccMoTpeHa 4acToTa Puna 
B JIMHeHHOŃ ONHO — A MHOTOMEpHOJi CHCTEME, OIMACAHHOi ypaBHeHHEM COCTOXHHH. B KOHNE 
pa60Tki upe4CTaBIeHKI pe3YJIbTATKI KOMIIBIOTEDHKIX HMATAIHA. 
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